s=10;b=8/3; r=28;h=0.003; a=28
x=0;y=1; z=0; f=zeros(1,3); f=[x y z];
k=1;
for t=0:h:100
k1_x=a*(y(k)-x(k));
k2_x=a*(y(k)-x(k))+(k1_x*h*0.5);
k3_x=a*(y(k)-x(k))+(k2_x*h*0.5);
k4_x=a*(y(k)-x(k))+(k3_x*h);
x(k+1)=x(k)+(k1_x+2*k2_x+2*k3_x+k4_x)*h/6;
k1_y=x(k)*(r-z(k))-y(k);
k2_y=(x(k)*(r-z(k))-y(k))+k1_y*h*0.5;
k3_y=(x(k)*(r-z(k))-y(k))+k2_y*h*0.5;
k4_y=(x(k)*(r-z(k))-y(k))+k3_y*h*0.5;
y(k+1)=y(k)+(k1_y+(2*k2_y)+(2*k3_y)+(k4_y))*h/6;
k1_z=x(k)*y(k)-b*z(k);
k2_z=(x(k)*y(k)-b*z(k))+k1_z*h*0.5;
k3_z=(x(k)*y(k)-b*z(k))+k2_z*h*0.5;
k4_z=(x(k)*y(k)-b*z(k))+k3_z*h*0.5;
z(k+1)=z(k)+(k1_z+(2*k2_z)+(2*k3_z)+(k4_z))*h/6;
k=k+1;
end
plot3(x,y,z); grid on; xlabel('x'); ylabel('y');
zlabel('z'); title('Lorenz attractor');
Best Answer