In[]:=
tet=10;​​r=28;​​b=0.1;​​e=(r-1)^(-1/2);​​a=e*b*tet^(-1/2);​​lmb=e*(b*tet+1)^(-1/2);​​Print["e=",e]​​Print["a=",a]​​Print["lmb=",lmb]​​xk=2;​​yk=1;​​zk=1;​​h=0.1;​​fx[t_,x_,y_,z_]=y;​​fy[t_,x_,y_,z_]=x-lmb*y-x*z;​​fz[t_,x_,y_,z_]=-a*z+x^2;
e=
1
3
3
a=0.00608581
lmb=0.136083
In[]:=
solution=NDSolve[{x'[t]==y[t],y'[t]==x[t]-lmb*y[t]-x[t]*z[t],z'[t]==-a*z[t]+x[t]^2,x[0]==2,y[0]==1,z[0]==1},{x,y,z},{t,100}];
In[]:=
Plot[{Evaluate[x[t]/.solution],Evaluate[y[t]/.solution],​​Evaluate[z[t]/.solution]},{t,0,100},PlotLegendsLineLegend[{x,y,z}]]
Out[]=
20
40
60
80
100
-5
5
10
15
x
y
z
In[]:=
DoPrint["Step h=",h];​​res=NSolve​​k1x==h*fxt+
1
2
-
3
6
h,xk+
1
4
k1x+
1
4
-
3
6
k2x,yk+
1
4
k1y+
1
4
-
3
6
k2y,zk+
1
4
k1z+
1
4
-
3
6
k2z,​​ k2x==h*fxt+
1
2
+
3
6
h,xk+
1
4
+
3
6
k1x+
1
4
k2x,yk+
1
4
+
3
6
k1y+
1
4
k2y,zk+
1
4
+
3
6
k1z+
1
4
k2z,​​k1y==h*fyt+
1
2
-
3
6
h,xk+
1
4
k1x+
1
4
-
3
6
k2x,yk+
1
4
k1y+
1
4
-
3
6
k2y,zk+
1
4
k1z+
1
4
-
3
6
k2z,​​ k2y==h*fyt+
1
2
+
3
6
h,xk+
1
4
+
3
6
k1x+
1
4
k2x,yk+
1
4
+
3
6
k1y+
1
4
k2y,zk+
1
4
+
3
6
k1z+
1
4
k2z,​​k1z==h*fzt+
1
2
-
3
6
h,xk+
1
4
k1x+
1
4
-
3
6
k2x,yk+
1
4
k1y+
1
4
-
3
6
k2y,zk+
1
4
k1z+
1
4
-
3
6
k2z,​​ k2z==h*fzt+
1
2
+
3
6
h,xk+
1
4
+
3
6
k1x+
1
4
k2x,yk+
1
4
+
3
6
k1y+
1
4
k2y,zk+
1
4
+
3
6
k1z+
1
4
k2z,​​{k1x,k2x,k1y,k2y,k1z,k2z},Reals;​​xk1=xk+
h
2
(res[[1,1,2]]+res[[1,2,2]]);​​yk1=yk+
h
2
(res[[1,3,2]]+res[[1,4,2]]);​​zk1=zk+
h
2
(res[[1,5,2]]+res[[1,6,2]]);​​Print["xk: ",xk1];​​Print["yk: ",yk1];​​Print["zk: ",zk1];​​xk=xk1;​​yk=yk1;​​zk=zk1;​​solution=NDSolve[{x'[t]==y[t],y'[t]==x[t]-lmb*y[t]-x[t]*z[t],z'[t]==-a*z[t]+x[t]^2,x[0]==2,y[0]==1,z[0]==1},{x,y,z},{t,100}];​​realx=x[h]/.solution;​​realy=y[h]/.solution;​​realz=z[h]/.solution;​​Print["real x=",realx];​​Print["real y=",realy];​​Print["real z=",realz];​​Print["fault x=",Abs[realx-xk1]];​​Print["fault y=",Abs[realy-yk1]];​​Print["fault z=",Abs[realz-zk1]];​​,{h,0,100,0.1}
Step h=0.
xk: 2.
yk: 1.
zk: 1.
real x={2.}
real y={1.}
real z={1.}
fault x={0.}
fault y={0.}
fault z={0.}
Step h=0.1
xk: 2.00979
yk: 0.994404
zk: 1.04194
real x={2.09793}
real y={0.94404}
real z={1.41936}
fault x={0.088136}
fault y={0.0503636}
fault z={0.377426}
Step h=0.2
xk: 2.04633
yk: 0.949183
zk: 1.21891
real x={2.18565}
real y={0.793453}
real z={1.87769}
fault x={0.139327}
fault y={0.15573}
fault z={0.658776}
Step h=0.3
xk: 2.11084
yk: 0.762287
zk: 1.6418
real x={2.25302}
real y={0.534905}
real z={2.36994}
fault x={0.142172}
fault y={0.227382}
fault z={0.72814}
Step h=0.4
xk: 2.14177
yk: 0.196248
zk: 2.41141
real x={2.28875}
real y={0.160053}
real z={2.88546}
fault x={0.146974}
fault y={0.0361949}
fault z={0.474046}
Step h=0.5
xk: 1.91527
yk: -1.0156
zk: 3.42668
real x={2.28116}
real y={-0.330294}
real z={3.40752}
fault x={0.365889}