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,
]; ]; ];
end
Lyapunov1 n=1:20001;
plot(n,Lya1,'k',n,Lya2,'k',n,Lya3,'k') %grid on
axis([0,30001,-0.8,0.5])
title('Lyapunov exponents of Warship') xlabel('n'),ylabel('LCE')