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]

哪位高手帮我瞧瞧，非常感谢！