求最大李雅普诺夫指数地matlab程序

b*cos(gama*(z(n)+h/2)).*sin(x(n)+k1/2)+c);  k3=h*(y(n)+m2/2);  m3=h*(-sin(x(n)+k2/2)-a*(y(n)+m2/2)+b*cos(gama*(z(n)+h/2)).*sin(x(n)+k2/2)+c);  k4=h*(y(n)+m3);  m4=h*(-sin(x(n)+k3)-a*(y(n)+m3)+b*cos(gama*(z(n)+h)).*sin(x(n)+k3)+c);  x(n+1)=x(n)+(k1+2*k2+2*k3+k4)/6;  y(n+1)=y(n)+(m1+2*m2+2*m3+m4)/6;  z(n+1)=n*h;  J=[0  1  0;      b*cos(gama*z(n+1))*cos(x(n+1))-cos(x(n+1))  -a  -b*gama*sin(gama*z(n+1))*sin(x(n+1));      0  0  0];  J=eye(3)+h*J;    B=J*V*S;    [V,S,U]=svd(B);  a_max=max(diag(S));  S=(1/a_max)*S;  b1=b1+log(a_max);  Lyapunov1=(log(diag(S))+b1)/(n*h);  Lya1=[Lya1,Lyapunov1(1,];  Lya2=[Lya2,Lyapunov1(2,];  Lya3=[Lya3,Lyapunov1(3,];endLyapunov1n=1:20001;plot(n,Lya1,'k',n,Lya2,'k',n,Lya3,'k')%gridonaxis([0,30001,-0.8,0.5])title('LyapunovexponentsofWarship')xlabel('n'),ylabel('LCE')