Методи на Рунге-Кута от II-ри ред
задача 2 а) от стр.35

Рунге-Кута от II-ри ред -  формула (1,1)

x0=1.; xend=2;n=4;y0=0;
h=((xend-x0)/n);
f[x_,y_]:=(y/x)+1;
yt[x_]:=x Log[x]; (* въвеждаме точното решение *)
x=x0;
y=y0;
points = {{x,y}};
Print["i = ",0, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y],  " \!\(\*SubscriptBox[\(y\), \(точно\)]\)(x) = ",yt[x], " грешка = ", Abs[y-yt[x]] ]
For[i=1,i≤ n,i++,
(* тяло на цикъла *)
k1 = h*f[x,y];
k2 = h*f[x+h,y+k1];
y=y+(1/2)(k1+k2);
x=x+h;
AppendTo[points,{x,y}];
Print["i = ",i, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y]," k1 = ",k1," k2 = ",k2,  " \!\(\*SubscriptBox[\(y\), \(точно\)]\)(x) = ",yt[x], " грешка = ", Abs[y-yt[x]] ]
]
grfyt=Plot[yt[x],{x,x0,xend}, PlotStyle->Red];
grfp=ListPlot[points, PlotStyle->{Black, PointSize[0.02]}];
grfmO=Show[grfyt,grfp]

Рунге-Кута от II-ри ред -  формула (1/2,1/2)

x0=1.; xend=2;n=4;y0=0;
h=((xend-x0)/n);
f[x_,y_]:=(y/x)+1;
yt[x_]:=x Log[x]; (* въвеждаме точното решение *)
x=x0;
y=y0;
points = {{x,y}};
Print["i = ",0, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y],  " \!\(\*SubscriptBox[\(y\), \(точно\)]\)(x) = ",yt[x], " грешка = ", Abs[y-yt[x]] ]
For[i=1,i≤ n,i++,
(* тяло на цикъла *)
k1 = h*f[x,y];
k2 = h*f[x+(h/2),y+(k1/2)];
y=y+k2;
x=x+h;
AppendTo[points,{x,y}];
Print["i = ",i, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y]," k1 = ",k1," k2 = ",k2,  " \!\(\*SubscriptBox[\(y\), \(точно\)]\)(x) = ",yt[x], " грешка = ", Abs[y-yt[x]] ]
]
grfyt=Plot[yt[x],{x,x0,xend}, PlotStyle->Red];
grfp=ListPlot[points, PlotStyle->{Black, PointSize[0.02]}];
grfmO=Show[grfyt,grfp]

Рунге-Кута от II-ри ред -  формула (2/3,2/3)

x0=1.; xend=2;n=4;y0=0;
h=((xend-x0)/n);
f[x_,y_]:=(y/x)+1;
yt[x_]:=x Log[x]; (* въвеждаме точното решение *)
x=x0;
y=y0;
points = {{x,y}};
Print["i = ",0, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y],  " \!\(\*SubscriptBox[\(y\), \(точно\)]\)(x) = ",yt[x], " грешка = ", Abs[y-yt[x]] ]
For[i=1,i≤ n,i++,
(* тяло на цикъла *)
k1 = h*f[x,y];
k2 = h*f[x+(2/3)h,y+(2/3)k1];
y=y+(1/4)k1 + (3/4)k2;
x=x+h;
AppendTo[points,{x,y}];
Print["i = ",i, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y]," k1 = ",k1," k2 = ",k2,  " \!\(\*SubscriptBox[\(y\), \(точно\)]\)(x) = ",yt[x], " грешка = ", Abs[y-yt[x]] ]
]
grfyt=Plot[yt[x],{x,x0,xend}, PlotStyle->Red];
grfp=ListPlot[points, PlotStyle->{Black, PointSize[0.02]}];
grfmO=Show[grfyt,grfp]

Рунге-Кута от IV-ри ред

x0=1.; xend=2;n=4;y0=0;
h=((xend-x0)/n);
f[x_,y_]:=(y/x)+1;
yt[x_]:=x Log[x]; (* въвеждаме точното решение *)
x=x0;
y=y0;
points = {{x,y}};
Print["i = ",0, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y],  " \!\(\*SubscriptBox[\(y\), \(точно\)]\)(x) = ",yt[x], " грешка = ", Abs[y-yt[x]] ]
For[i=1,i≤ n,i++,
(* тяло на цикъла *)
k1 = h*f[x,y];
k2 = h*f[x+(1/2)h,y+(1/2)k1];
k3 = h*f[x+(1/2)h,y+(1/2)k2];
k4 = h*f[x+h,y+k3];
y=y+(1/6)(k1 + 2k2+ 2k3+k4);
x=x+h;
AppendTo[points,{x,y}];
Print["i = ",i, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y]," k1 = ",k1," k2 = ",k2,  " \!\(\*SubscriptBox[\(y\), \(точно\)]\)(x) = ",yt[x], " грешка = ", Abs[y-yt[x]] ]
]
grfyt=Plot[yt[x],{x,x0,xend}, PlotStyle->Red];
grfp=ListPlot[points, PlotStyle->{Black, PointSize[0.02]}];
grfmO=Show[grfyt,grfp]