s=10; %sigma b=8/3; %beta r=28; %rho h=0.003; %time step size x=0; %initialize x y=1; z=0; f=zeros(1,3); f=[x y z]; for t=0:h:100 k1_x=a*(y-x); k2_x=a*(y-x)+(k1_x*h*0.5); k3_x=a*(y-x)+(k2_x*h*0.5); k4_x=a*(y-x)+(k3_x*h); x_n=x+(k1_x+(2*k2_x)+(2*k3_x)+(k4_x))*h/6; %x_n is xn+1
k1_y=x*(r-z)-y; k2_y=(x*(r-z)-y)+k1_y*h*0.5; k3_y=(x*(r-z)-y)+k2_y*h*0.5; k4_y=(x*(r-z)-y)+k3_y*h*0.5; y_n=y+(k1_y+(2*k2_y)+(2*k3_y)+(k4_y))*h/6; %y_n is yn+1, the next value of y
k1_z=x*y-b*z; k2_z=(x*y-b*z)+k1_z*h*0.5; k3_z=(x*y-b*z)+k2_z*h*0.5; k4_z=(x*y-b*z)+k3_z*h*0.5; z_n=z+(k1_z+(2*k2_z)+(2*k3_z)+(k4_z))*h/6; %z_n is zn+1 the next value of z
x=x_n; y=y_n; z=z_n;
end; plot3(x,y,z); grid on; xlabel('x'); ylabel('y'); zlabel('z'); title('Lorenz attractor');
The plot comes back empty.
Best Answer