Програмен код, който реализира метода на Ойлер за решаване на ОДУ:

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;
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++,
(* тяло на цикъла *)
y=y+h*f[x,y];
x=x+h;
Print["i = ",i, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y],  " \!\(\*SubscriptBox[\(y\), \(точно\)]\)(x) = ",yt[x], " грешка = ", Abs[y-yt[x]] ]
]

x0=1.; xend=2;n=10;y0=0;
h=((xend-x0)/n);
f[x_,y_]:=(y/x)+1;
yt[x_]:=x Log[x]; (* въвеждаме точното решение *)
x=x0;
y=y0;
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++,
(* тяло на цикъла *)
y=y+h*f[x,y];
x=x+h;
Print["i = ",i, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y],  " \!\(\*SubscriptBox[\(y\), \(точно\)]\)(x) = ",yt[x], " грешка = ", Abs[y-yt[x]] ]
]

x0=1.; xend=2;n=100;y0=0;
h=((xend-x0)/n);
f[x_,y_]:=(y/x)+1;
yt[x_]:=x Log[x]; (* въвеждаме точното решение *)
x=x0;
y=y0;
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++,
(* тяло на цикъла *)
y=y+h*f[x,y];
x=x+h;
Print["i = ",i, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y],  " \!\(\*SubscriptBox[\(y\), \(точно\)]\)(x) = ",yt[x], " грешка = ", Abs[y-yt[x]] ]
]

визуализация на получените резултати във вид на точки и точното решение:

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++,
(* тяло на цикъла *)
y=y+h*f[x,y];
x=x+h;
AppendTo[points,{x,y}];
(*Print["i = ",i, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y],  " \!\(\*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]}];
Show[grfyt,grfp]

x0=1.; xend=2;n=10;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++,
(* тяло на цикъла *)
y=y+h*f[x,y];
x=x+h;
AppendTo[points,{x,y}];
(*Print["i = ",i, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y],  " \!\(\*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]}];
Show[grfyt,grfp]

x0=1.; xend=2;n=100;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++,
(* тяло на цикъла *)
y=y+h*f[x,y];
x=x+h;
AppendTo[points,{x,y}];
(*Print["i = ",i, " x = ",x, " y = ",y, " f(x,y) = ",f[x,y],  " \!\(\*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.01]}];
Show[grfyt,grfp]