洛伦兹方程是一组典型的非线性方程
for i=1:1:3000
abx=sym('10*(y0-x0)');
aby=sym('x0*(28-z0)-y0');
abz=sym('x0*y0-8/3*z0');
x0=x1;y0=y1;z0=z1;
kfx=eval_r(abx);kfy=eval_r(aby);kfz=eval_r(abz);
x0=x1+1/2*h*kfx;y0=y1+1/2*h*kfy;z0=z1+1/2*h*kfz;
ksx=eval_r(abx);ksy=eval_r(aby);ksz=eval_r(abz);
x0=x1+1/2*h*ksx;y0=y1+1/2*h*ksy;z0=z1+1/2*h*ksz;
ktx=eval_r(abx);kty=eval_r(aby);ktz=eval_r(abz);
x0=x1+1/2*h*ktx;y0=y1+1/2*h*kty;z0=z1+1/2*h*ktz;
kfhx=eval_r(abx);kfhy=eval_r(aby);kfhz=eval_r(abz);
x2=x1+h/6*(kfx+2*ksx+2*ktx+kfhx);
y2=y1+h/6*(kfy+2*ksy+2*kty+kfhy);
z2=z1+h/6*(kfz+2*ksz+2*ktz+kfhz);
a(i)=x2;b(i)=y2;c(i)=z2;
x1=x2;y1=y2;z1=z2;
end
plot3(a,b,c)