求助：怎么在M文件里用求出非线性函数的数值解

Cavatina

源程序如下：
function [epsilon,miu]=retrieve(E_0,E_r,l,w,theta)
%retrieve the equivalent paremeters of matamaterials
%E_0 is the input field at the original point with zero phase
%E_r is the reflectd field with certain phase
%l is the length of the material
%epsilon=ε；miu=μ;ommiga=w;gamma=γ_12 is the ratio of reflection at the 1-2border
%start calculate
%subtitude the angle of the complex exponetial function:F=-2*l*k_z;
[F,gamma]=solve('real(E_r/E_0)-((1-cos(F))*(1+gamma^2))/(gamma^4+2*cos(F)*gamma^2+1)=0','imga(E_r/E_0)-(sin(F)*(gamma^2-1))/(gamma^4+2*cos(F)*gamma^2+1)=0','F,gamma');
subs([F,gamma]);
%Y_1 is the intrinsic admittance of the air
k_z=F/(l*(-2));
Y_1=(8.85*10^(-5)/(4*3.1416))^(1/2);
miu=solve('gamma=(w*miu*Y_1-k_z)/(w*miu*Y_1+k_z)','miu');
%using the dispersion formula to calculate epsilon
k=sqrt((k_z).^2+(w*sin(theta)).^2*8.85*4*3.1416*10^(-19));
epsilon=k.^2/(w.^2*miu);
end
上面是之前写的一个计算电磁材料等效系数的函数文件，里面大部分内容其实是备注解释，有用的东西就是solve里面的函数，主要遇到了以下问题：
1.所求函数为复函数，其中F和gamma都是实数，但是因为方程包含实部和虚部，所以我用real和imga分开了，不过貌似求不来
2.一旦在M文件里用solve函数，里面包含的变量在系统看来是没用过的（他会提示我未使用变量）
3.有人建议使用subs赋值求解，不过那样相当于先求出解析解，我试了一下，好像方程太复杂程序直接崩溃了

  我寒假自学的matlab，这么个简单的问题困扰我一个月了，我觉得浪费了很多时间，希望高手一次性给个痛快！！（方便的话最好给个改过的例子。。。感觉很多方法懂和实际写还是有差距）