CITE THIS NOTEBOOK: Introduction to numerical methods: discretization of the continuous by Athanasios Paraskevopoulos. Wolfram Community JAN 27 2023.
I am going to solve some exercises with the Classic Euler’s Method and the Improved Euler’s method for first-order equations.
The following exercise set is from the book Mathematica and Applications for Mathematicians, Physicists, and Engineers by Stephanos Trachanas.

FIRST-ORDER EQUATIONS

In[]:=
Clear["Global`*"];
1
.
Classic Euler' s method :
In[]:=
P[h_,N_]:=(u[0]=1;Do[u[n+1]=u[n]+h*f[n*h,u[n]],{n,0,N}])
2
.
Improved Euler' s Method:
In[]:=
Q[h_,N_]:=(u[0]=1;Do[u[n+1]=u[n]+h*f[n*h+(h/2),u[n]+(h/2)*f[n*h,u[n]]],{n,0,N}])
3
.
Both Methods in a single command :
In[]:=
P[q_,h_,N_]:=(u[0]=1;Do[u[n+1]=u[n]+h*f[n*h+(h*q/2),u[n]+(h*q/2)*f[n*h,u[n]]],{n,0,N}])
When q=1 we use Improved Euler’ s Method and when q=0 we use Classic Euler’ s method.
h : step of partitioning of the interval
N : the number of points

Example 1

Solve y' (x) = 3 Exp[-x] - y (x), y (0) = 1. We will study this function on the interval [0 , 5]
Clear["Global`*"];
In[]:=
f[x_,y_]=3Exp[-x]-y
Out[]=
3
-x

-y
Classic Euler result :
In[]:=
P[0,0.1,50];
In[]:=
u[50]
Out[]=
0.103399
Improved Euler result :
In[]:=
P[1,0.1,50]
In[]:=
u[50]
Out[]=
0.107909
Exact solution :
In[]:=
y[x_]=(3x+1)Exp[-x]
Out[]=
-x

(1+3x)
In[]:=
y[5.]
Out[]=
0.107807
In[]:=
Y=Interpolation[Table[{0.1n,u[n]},{n,0,50}]]
Out[]=
InterpolatingFunction
Domain: {{0.,5.}}
Output: scalar

In[]:=
Plot[Y[x],{x,0,5}]
Out[]=
1
2
3
4
5
0.5
1.0
1.5

Example 2

Solve y’(x)=Exp[-x*y] , y(0)=1. We will study this function on the interval [0 , 3] with step h=0.01
Clear["Global`*"];
In[]:=
f[x_,y_]=Exp[-x*y]
Out[]=
-xy

Classic Euler' s result :
In[]:=
P[0,0.01,300]
In[]:=
u[300]
Out[]=
1.66427
Improved Euler' s result :
In[]:=
P[1,0.01,300]
In[]:=
u[300]
Out[]=
1.66045
In[]:=
Y=Interpolation[Table[{0.01n,u[n]},{n,0,300}]]
Out[]=
InterpolatingFunction
Domain: {{0.,3.}}
Output: scalar

In[]:=
Plot[Y[x],{x,0,3}]
Out[]=
0.5
1.0
1.5
2.0
2.5
3.0
1.1
1.2
1.3
1.4
1.5
1.6
Comparison with the NDSolve result :
In[]:=
Clear[y]
In[]:=
s=NDSolve[{y'[x]Exp[-x*y[x]],y[0]1},y,{x,0,3}];
In[]:=
y[3]/.s
Out[]=
{1.66046}
What do you notice?

Example 3

In[]:=
Clear["Global`*"];
Solve y' (x) = Sin[x*y] , y (0) = 1. We will study this function on the interval [0 , 10] with step h = 0.01
In[]:=
f[x_,y_]=Sin[x*y]
Out[]=
Sin[xy]
Classic Euler' s result :
In[]:=
P[0,0.01,1000]
In[]:=
u[1000]
Out[]=
0.317397
Improved Euler' s result :
In[]:=
P[1,0.01,1000]
In[]:=
u[1000]
Out[]=
0.317401
In[]:=
Y=Interpolation[Table[{0.01n,u[n]},{n,0,1000}]]
Out[]=
InterpolatingFunction
Domain: {{0.,10.}}
Output: scalar

In[]:=
Plot[Y[x],{x,0,10}]
Out[]=
2
4
6
8
10
0.5
1.0
1.5
2.0
In[]:=
Table[Y[n],{n,1,10}]
Comparison with the NDSolve result :

Example 4

Solve y' (x) = x^2 - y^2 , y (0) = 1. We will study this function on the interval [0 , 3] with step h = 0.01
Classic Euler' s result :
Improved Euler' s result :
Comparison with NDSolve result :