matlab逆求贝塞尔函数变量值,MATLAB怎么求解有贝塞尔函数的问题，求高手帮帮忙，谢谢...
[code]%%声波声速、声衰减计算公式程序%%%%%
syms kb u w n p q m kesi x
coef={kb,u,w,n,p,q,m,kesi};
%u为骨架剪切模量,kb为骨架体变模量,p=(1-n)ps固相多孔隙介质的等效密度,q=npf流体的等效密度,w圆
频率,kr颗粒体变模量,kf流体体变模量,nta粘滞系数,kao渗透率,alfa弯曲度,a孔隙尺寸,n孔隙度%%；
a=5e-5;
kr=3.6e10;
kf=2.25e9;
nta=0.001;
kao=1.0e-10;
alfa=1.25;
D=kr*(1+(kr/kf-1)*n);
H=(kr-kb)^2/(D-kb)+kb+4*u/3;
C=kr*(kr-kb)/(D-kb);
M=kr^2/(D-kb);
m=1.25*1023/n;
for f=1:1200
tt(f)=f;
w=2*pi*f;
kesi=a*(w*1023/nta).^1/2;
T=besselj(1,kesi)./besselj(0,kesi);
F=kesi.*T/4/(1-2*T/(i.*kesi));
A=C^2-H*M;
B=m*H*w^2+1999.2*w^2*M-i*H*w*F*nta/kao-2*C*1023*w^2;
C=-1999.2*m*w^4+i*1999.2*F*nta*w^3/kao+1023^2*w^4;
equa=A*x^4+B*x^2+C;
g=solve(equa,'x');
n=0.5;
lr=subs(g,coef,{4.4e7-i*2.0e6,2.6e7-i*1.25e6,2*pi*f,n,(1-n)*2650,n*1023,1.25*1023/n,5.0e-5*
(2*pi*f*1023/0.001)^1/2});
E=real(lr);
vp=[1 1 1 1]'./E;
li=subs(-w*g,coef,{4.4e7-i*2.0e6,2.6e7-i*1.25e6,2*pi*f,n,(1-n)*2650,n*1023,1.25*1023/n,5.0e
-5*(2*pi*f*1023/0.001)^1/2});
F1=imag(li);
ap=F1;
v=[vp(1) 0];
a=[ap(1) 0];
for j=2:4
if vp(j)>vp(1)
v(1)=vp(j);
a(1)=ap(j);
end
end
v1(f)=v(1);
a1(f)=a(1);
end
subplot(1,2,1);plot(v1,'-','linewidth',2);xlabel('f/Hz');ylabel('Vp/m/s');title('沉积物声速-
频率计算曲线');
subplot(1,2,2);plot(a1,'-','linewidth',2);xlabel('f/Hz');ylabel('ap/dB/m');title('沉积物声衰
减-频率计算曲线');
Function 'gt' is not implemented for MuPAD symbolic objects.[
不知道程序怎么弄了，有两个问题1.引入了贝塞尔函数，2.有复数解情况/code]
哪位高手帮我瞧瞧，非常感谢！