Euler’s method is one of the simplest numerical methods for solving initial value problems. In this section, we discuss the theory and implementation of Euler’s method in matlab.

Leonhard Euler was born in 1707, Basel, Switzerland and passed away in 1783, Saint Petersburg, Russia. In 1738, he became almost blind in his right eye. Euler was one of the most eminent mathematicians of the 18th century, and is held to be one of the greatest in history. He is also widely considered to be the most prolific mathematician of all time. He spent most of his adult life in Saint Petersburg, Russia, except about 20 years in Berlin, then the capital of Prussia.

In 1768, Leonhard Euler (St. Petersburg, Russia) introduced a numerical method that is now called the Euler method or the tangent line method for solving numerically the initial value problem:

\[ y' = f(x,y), \qquad y(x_0 ) = y_0 , \]
where f(x,y) is the instantaneous slope function, and \( (x_0 , y_0 ) \) is a prescribed point on the plane. Euler's method or rule is a very basic algorithm that could be used to generate a numerical solution to the initial value problem for first order differential equation. The solution that it produces will be returned to the user in the form of a list of points. The user can then do whatever one likes with this output, such as create a graph, or utilize the point estimates for other purposes. Euler's rule serves to illustrate the concepts involved in the advanced methods. It has limited usage because of the larger error that is accumulated as the process proceeds, so it requires a smaller step size. However, it is importatnt to study because the error analysis is easier to understand.

To start, we need mesh or grid points, that is, the set of discrete points for an independent variable at which we find approximate solutions. In other words, we will find approximate values of the unknown solution at these mesh points. For simplicity, we use uniform grid with fixed step length h; however, in practical applications, it is almost always not a constant. For convenience, we subdivide the interval of interest [𝑎,b] (where usually 𝑎 = x0 is the starting point) into m equal subintervals and select the mesh/grid points

\[ x_n = a + n\,h \qquad\mbox{for} \quad n=0,1,2,\ldots , m; \quad\mbox{where} \quad h = \frac{b-a}{m} . \]
The considered initial value problem is assumed to have a unique solution \( y= \phi (x) , \) on the interval of interest {a,b], and its approximations at the grid points will be denoted by yn, so we wish that \( y_n \approx \phi (x_n ) , \quad n=1,2, \ldots . \)

There are three main approaches (we do not discuss others in this section) to derive the Euler rule: either use finite difference approximation to the derivative or transfer the initial value problem to the Volterra integral equation and then truncate it, or apply Taylor series.

  1. We start with the finite difference method. It is known that numerical differentiation is an ill-posed problem that does not depend continuously on the data. If we approximate the derivative in the left-hand side of the differential equation y' = f(x,y) by the finite difference

    \[ y' (x_n ) \approx \frac{y_{n+1} - y_n}{h} \]
    on the small subinterval \( [x_{n+1} , x_n ] , \) we arrive at the Euler's rule when the slope function is evaluated at x = xn.
    \begin{equation*} y_{n+1} = y_n + (x_{n+1} - x_n ) f( x_n , y_n ) \qquad \mbox{or} \qquad y_{n+1} = y_n + h f_n , \end{equation*}
    where the following notations are used: \( h=x_{n+1} - x_n \) is the step length (which is assumed to be constant for simplicity), \( f_n = f( x_n , y_n ) \) is the value of slope function at mesh point, and yn denotes the approximate value of the actual solution y = ϕ(x) at the point \( x_n = x_0 + n\,h \) (\( n=1,2,\ldots \) ).
  2. Integrating the identity y'(x) = f(x,y(x)) between xn and xn+1, we obtain
       
    curve = Plot[x - Exp[x - 0.05] + 1.5, {x, -1.0, 0.7}, PlotStyle -> Thick, Axes -> False, Epilog -> {Blue, PointSize@Large, Point[{{-0.5, 0.42}, {0.25, 0.53}}]}, PlotRange -> {{-1.6, 0.6}, {-0.22, 0.66}}];;
    line = ListLinePlot[{{-0.5, 0.42}, {0.25, 0.42}}, Filling -> Bottom, FillingStyle -> Opacity[0.6]];
    curve2 = Plot[x - Exp[x - 0.05] + 1.5, {x, -0.5, 0.25}, PlotStyle -> Thick, Axes -> False, PlotRange -> {{-1.6, 0.6}, {0.0, 0.66}}, Filling -> Bottom, FillingStyle -> Yellow];;
    ar1 = Graphics[{Black, Arrowheads[0.07], Arrow[{{-1.0, 0.0}, {0.6, 0.0}}]}];;
    t1 = Graphics[ Text[Style["f(t, y(t))", FontSize -> 12], {-0.7, 0.5}]]; xn = Graphics[ Text[Style[Subscript["x", "n"], FontSize -> 12, FontColor -> Black, Bold], {-0.48, -0.05}]];;
    xn1 = Graphics[ Text[Style[Subscript["x", "n+1"], FontSize -> 12, FontColor -> Black, Bold], {0.24, -0.05}]];;
    Show[curve, curve2, line, ar1, t1, xn, xn1]
    \[ y(x_{n+1}) - y(x_n ) = \int_{x_n}^{x_{n+1}} f(t, y(t))\,{\text d} t . \]
    Such integration transfers the initial value problem for unbounded derivative operator into the integral equation, avoiding a discontinuous operator. Approximating the integral by the left rectangular (very crude) rule for numerical integration (length of interval times value of integrand at left end point) and identifying \( y(x_n ) \) with yn, we again obtain the Euler formula.

  3. We may assume the possibility of expanding the solution in a Taylor series around the point xn:
    \[ y(x_{n+1}) = y(x_n +h) = y(x_{n}) + h\, f(x_n, y(x_n )) + \frac{1}{2}\, h^2 \,y'' (x_n ) + \cdots . \]
    The Euler formula is the result of truncating this series after the linear term in h:    yn+1 = yn + h f(xn, yn).
curve = Plot[x - Exp[x - 0.05] + 1.5, {x, -1.0, 0.7}, PlotStyle -> Thick, Axes -> False,
Epilog -> {Blue, PointSize@Large, Point[{{-0.4, 0.46}, {-0.4, 0.5}, {0, 0.63}, {0, 0.55},
{0.4, 0.48}, {0.4, 0.61}}]}, PlotRange -> {{-1.6, 0.6}, {-0.22, 0.66}}]
line1 = Graphics[Line[{{-0.8, 0.27}, {-0.4, 0.5}}]];
line2 = Graphics[{Dashed, Line[{{-0.8, 0.27}, {-0.8, -0.16}}]}];
line2a = Graphics[{Dashed, Line[{{-0.8, 0.27}, {-1.5, 0.27}}]}];
line3 = Graphics[{Dashed, Line[{{-0.4, 0.46}, {-0.4, -0.16}}]}];
line3a = Graphics[{Dashed, Line[{{-0.4, 0.46}, {-1.5, 0.46}}]}];
line4 = Graphics[{Dashed, Line[{{0.0, 0.55}, {0.0, -0.16}}]}];
line4a = Graphics[{Dashed, Line[{{0.0, 0.55}, {-1.5, 0.55}}]}];
line5 = Graphics[{Dashed, Line[{{0.4, 0.48}, {0.4, -0.16}}]}];
ar1 = Graphics[{Black, Arrowheads[0.07], Arrow[{{-1.5, -0.22}, {-1.5, 0.65}}]}];
ar2 = Graphics[{Black, Arrowheads[0.07], Arrow[{{-1.5, -0.16}, {0.6, -0.16}}]}];
aa1 = Graphics[{Arrowheads[{-0.04, 0.04}], Arrow[{{-0.8, 0}, {-0.4, 0}}]}];
aa2 = Graphics[{Arrowheads[{-0.04, 0.04}], Arrow[{{-0.4, 0}, {-0.0, 0}}]}];
aa3 = Graphics[{Arrowheads[{-0.04, 0.04}], Arrow[{{0.0, 0}, {0.4, 0}}]}];
t1 = Graphics[Text[Style["h", FontSize -> 12], {-0.6, 0.05}]];
t2 = Graphics[Text[Style["h", FontSize -> 12], {-0.6, 0.05}]];
t3 = Graphics[Text[Style["h", FontSize -> 12], {0.2, 0.05}]];
x0 = Graphics[ Text[Style[Subscript["x", "0"], FontSize -> 12, FontColor -> Black, Bold], {-0.78, -0.19}]];
x1 = Graphics[ Text[Style[Subscript["x", "1"], FontSize -> 12, FontColor -> Black, Bold], {-0.38, -0.19}]];
x2 = Graphics[ Text[Style[Subscript["x", "2"], FontSize -> 12, FontColor -> Black, Bold], {0.02, -0.19}]];
x3 = Graphics[ Text[Style[Subscript["x", "3"], FontSize -> 12, FontColor -> Black, Bold], {0.42, -0.19}]];
x = Graphics[ Text[Style["x", FontSize -> 12, FontColor -> Black, Bold], {0.5, -0.09}]];
y = Graphics[ Text[Style["y", FontSize -> 12, FontColor -> Black, Bold], {-1.4, 0.6}]];
line6 = Graphics[Line[{{-0.4, 0.5}, {0.0, 0.63}}]];
y0 = Graphics[ Text[Style[Subscript["y", "0"], FontSize -> 12, FontColor -> Black, Bold], {-0.7, 0.26}]];
y1 = Graphics[ Text[Style[Subscript["y", "1"], FontSize -> 12, FontColor -> Black, Bold], {-0.47, 0.5}]];
y2 = Graphics[ Text[Style[Subscript["y", "2"], FontSize -> 12, FontColor -> Black, Bold], {0.09, 0.63}]];
y3 = Graphics[ Text[Style[Subscript["y", "3"], FontSize -> 12, FontColor -> Black, Bold], {0.35, 0.58}]];
tr = Graphics[ Text[Style["true value", FontSize -> 12, FontColor -> Black, Bold], {-1.2, 0.42}]];;
Show[curve, line1, line2, line3, line2a, line3a, line4, line4a , \ line5, ar1,
ar2, aa1, aa2, aa3, t1, t2, t3, x0, x1, x2, x3, x, y, \ line6, y0, y1, y2, y3,tr]
Each of these interpretations points the way to a class of generalizations of Euler's method that are discussed later. It is interesting to note that the generalization indicated by (i) which seems to be the most straightforward, has proved to be the least fruitful of the three.

In order to implement the Euler method, we need to follow the Euler algorithm:

  1. Read in the slope function and the initial values of all of the variables.
  2. Initialize the solution list to contain the initial condition point and define the step size or the number of steps.
  3. Determine next pair (xn+1, yn+1) according to the Euler rule:    yn+1 = yn + h f(xn, yn). Therefore, every next Euler point yn+1 is determined according to the tangent line starting from the previous grid point (xn,yn) with slope f(xn,yn); recall that the equation of the tangent line is y = yn + f(xn,yn) (x - xn).
\[ y_{n+1} = y_n + h\, f(x_{n}, y_{n}) , \qquad y_0 = y(0), \qquad n=0,1,2,\ldots . \]

When slope function is evaluated at the the right end point, we come to so called backward Euler rule:
\[ y_{n+1} = y_n + h\, f(x_{n+1}, y_{n+1}) , \qquad y_0 = y(0), \qquad n=0,1,2,\ldots . \]
We present the basic matlab code to implement Euler's rule
function  E=Euler(f,a,b,ya,m)
% Input: f is the handle of anonymous function
%        a and b are the left and right end points
%        ya is the initial condition y(a)
%        m is the number of steps
% Output: E=[x,y], where x is the vector of abscissas and y is the vector of ordinates
             
h=(b-a)/m;      % h is the step
y=zeros(m+1,1); % memory allocation, zero column vector
x=(a:h:b)';     % column vector of abscissas
y(1)=ya;
    for j=1:m
    y(j+1)=y(j)+h*f(x(j),y(j)); % Euler's formula
    end;
E=[x,y];
Its application is demonstrated in the following example.

Example: Consider the initial value problem \( y' = \left( x^2 - y^2 \right) \sin x, \qquad y(0) =1 . \) We present the direction field for the equation and result of calculations according to Euler's rule:

f=@ (x,y) (x^2-y^2)*sin(x) 
y=Euler(f,1,2,1,40); 
plot(y(:,1),y(:,2))
Another version:
function y=Euler1(n,t0,t1,x0)
% Input:        n is the number of steps
%               t0 and t1 are the left and right end points
%               x0 is the initial condition x(t0)
% Output:       y=[t,x]; (n+1) X 2 matrix
h=(t1-t0)/n;  % h is the step
t(n+1,1) = 0; % memory allocation for abscissas
x(n+1,1) = 0; % memory allocation for ordinates
t(1)=t0;
x(1)=x0;
for i=1:n
t(i+1)=t(1)+i*h;
x(i+1)=x(i)+h*Slope(t(i)); % Euler's formula
end;
y=[t,x]; % t and x are vectors of abscissas and ordinates respectively
Save this code as Euler1.m and prepare a file with slope function, say
\[ f(x) = \begin{cases} x, & \ \mbox{if } 0 \le x < 1, \\ 1, & \ \mbox{if } 1 \le x \le 2 , \\ 2- (x-1)^2 , & \ \mbox{if } x > 2. \end{cases} \]
Then we define this function in matlab as Slope

function f=Slope(x)
if x <=1
  f=x;
elseif 1 < x && x <=2
  f=1;
elseif x > 2
  f=2 - (x-1)^2;
end

Finally, we plot solution found with Euler's rule
N=60;
y=Euler1(N, 0, 3, 1); 
plot(h*(0:n),y(1:n+1),'linewidth',2)
title(['Euler Method using N=',num2str(N),' steps])
Now, on matlab prompt, you type Euler1(n,t0,t1,y0) and return, where n is the number of t-values, t0 and t1 are the left and right end points and y(t0) = y0 is the initial condition. You can accomplish this task without using function option, say

clear t % Clears old time steps and 
clear y % y values from previous runs 
a=0; % Initial time 
b=1; % Final time 
N=10; % Number of time steps 
y0=0; % Initial value y(a) 
h=(b-a)/N; % Time step 
t(1)=a; 
y(1)=y0; 
for n=1:N % For loop, sets next t,y values 
           t(n+1)=t(n)+h; 
           y(n+1)=y(n)+h*f(t(n),y(n)); % Calls the function f(t,y)=dy/dt 
end 
plot(t,y) title([’Euler Method using N=’,num2str(N),’ steps, by MYNAME’]) 

Another version:
function EulersMethod
%Eulers method
%y' = y/t - (y/t)^2 from t = [1 2] y(1) = 1 and h=0.1

a = 1;          %Interval we are using, from a to b.
b = 2;
y(1) = 1;       %Initial value of y, i.e. y = 1 for t = 0.
N = 10;         %Number of steps we will use
t(1) = a;       %inital time we will start at
h = (b - a) / N;    %Step size

F = @(y, t) y ./t - (y ./t).^2;  %model we are using

for i=1:N
    y(i+1) = y(i) + h*F(y(i), t(i));    %Implementation of Eulers method
    t(i+1) = a + i*h;
end

plot(t, y)

end

We start with finite difference method. If we approximate the derivative in the left hand-side by the finite difference

\[ y' (x_n ) \approx \frac{y_{n+1} - y_n}{h} \]
on the small subinterval \( [x_{n+1} , x_n ] , \) we arrive at the Euler's rule when right hand-side is evaluated at x = xn. When slope function is evaluated at the the right end point, we come to so called backward Euler rule:
\[ y_{n+1} = y_n + h\, f(x_{n+1}, y_{n+1}) , \qquad y_0 = y(0), \qquad n=0,1,2,\ldots . \]
However, this is not the only one approximation possible, one may consider the central difference:
\[ y' (x_n ) \approx \frac{y_{n+1} - y_{n-1}}{2h} , \]
which usually gives more accurate approximation. Using central difference, we get two-step approximation:
\[ y_{n+1} = y_{n-1} + 2h\, f\left( x_n , y_n \right) , \qquad n=1,2,\ldots , \]
which requires two starting points to solve this recurrence of second order. So one can use standard Euler rule to determine y1: \( y_1 = y_0 + h\, f\left( x_0 , y_0 \right) . \) However, this numerical algorithm may produce instability and we may observe numerical solution that is different from the true solution. Such solution is called "ghost solution." This henomenon is caused by roundoff errors.

To avoid the instability and to prevent the appearance of ghost solutions, many refine numerical procedures to integrate the ordinary differential equations have been proposed. In practical calculations, the central difference scheme is used when evaluations of slope function in intermediate points (as, for instance, required by Runge--Kutta methods) are not permitted. We mention two remedies to remove ghost solutions.

Consider the following mixed difference scheme parameterized by μ, 0 <μ<1, for the first order differential equation:
\[ \left( 1 - \mu \right) \frac{y_{n+1} - y_{n-1}}{2h} + \mu\,\frac{y_{n+1} - y_n}{h} = f\left( x_n , y_n \right) , \qquad n=1,2,\ldots . \]
If we set μ = 1, then we have Euler's rule. When μ =0, the scheme is nothing but the central difference scheme. This algorithm does not produce ghost solutions when μ > h/2.

Let φ be a real-valued function on \( \mathbb{R} \) that satisfies the property:

\[ \phi (h) = h + O\left( h^2 \right) \qquad\mbox{and} \qquad 0 < \phi(h) < 1 \quad\mbox{for all positive $h$}. \]
There exists a variety of functions φ that satisfy the above condition, e.g., \( \phi (h) = 1 - e^{-h} \quad \Longrightarrow \quad \phi (hq)/q = \left( 1- e^{-hq} \right) /q . \) Then the following nonstandard schemes are stable:
  • nonstandard explicit Euler method given by
    \[ \frac{y_{k+1} - y_k}{\phi (hq)/q} = f(y_k ), \qquad y_0 = y(0), \quad k=0,1,2,\ldots . \]
  • nonstandard implicit Euler method given by
    \[ \frac{y_{k+1} - y_k}{\phi (hq)/q} = f(y_{k+1} ), \qquad y_0 = y(0), \quad k=0,1,2,\ldots . \]

Example: Consider the initial value problem for the logistic equation

\[ \frac{{\text d}y}{{\text d}t} = y \left( 1-y \right) , \qquad y(0) = y_0 =0.5. \]
This problem has a unique solution
\[ y(t) = \frac{C\, e^t}{C\,e^t +1} , \qquad \mbox{where} \qquad C= \frac{y_0}{1-y_0} . \]
Using the following matlab code
function y=finite
y0=0.5;
h=.1;
n=500;
y(n)=0;
y1=y0+h*y0*(1-y0);
y(1)=y0;
y(2)=y1;
for i=3:n
    y(i)=y(i-2)+h*y(i-1)*(1-y(i-1));
end
plot(finite)
we observe instability

Another approach is based on transferring the given IVP \( y' = f(x,y) , \quad y(x_0 ) = y_0 \) to equivalent integral equation

\[ y(x_{n+1}) = y(x_n ) + \int_{x_n}^{x_{n+1}} f(s, y(s))\,{\text d} s , \qquad n=0,1,2,\ldots . \]
Application of left Riemann approximation, we come to the Euler's rule.

Example: To demonstrate the latter approach, we consider the initial value problem for the integro-differential equation

\[ \dot{y} = 2.3\, y -0.01\,y^2 -0.1\,y\, \int_0^t y(\tau )\,{\text d}\tau , \qquad y(0) =50. \]
Choosing the uniform grid \( t_k = k\,h , \quad k=0,1,2,\ldots , m; \) we integrate both sides to obtain the integral equation
\[ y(t ) = y(0) + 2.3\,\int_0^t y(s)\,{\text d}s -0.01\,\int_0^t y^2 (s) \,{\text d}s -0.1\,\int_0^t y(s)\,{\text d}s \int_0^s y(\tau )\,{\text d}\tau \]
for each mesh point t = tk. Since the double integral can be written as
\[ \int_0^t y(s)\,{\text d}s \int_0^s y(\tau )\,{\text d}\tau = \int_0^t \,{\text d}\tau \,\int_{\tau}^t {\text d}s \, y(s)\,y(\tau ) , \]
application of left rectangular rule yields for the first mesh point t1 = h:
\[ y(t_1 ) \approx y(0) + 2.3\,h\,y(0) -0.01\,h\,y^2 (0) -0.1 \, h\, \int_0^{t_1} {\text d}s \, y(s)\,y(0 ) \approx y(0) +2.3\,h\,y(0) -0.01\,h \,y(0)^2 - 0.1\,h^2 \,y(0)^2 . \]
So we choose the right hand-side as the approximate value y1.

For the general step in Euler's rule, we have

\[ y_{k+1} = y_k + 2.3\,h\,y_k -0.01\,h\,y^2_k -0.1 \, y_k \,h \int_0^{t_k} {\text d}s \, y(s)\qquad k=1,2,\ldots . \]
If the trapezoidal rule is used to approximate the integral, then this expression becomes
\[ y_{k+1} = y_k + 2.3\,h\,y_k -0.01\,h\,y^2_k -0.1 \, y_k \,h\, T_k (h) ,\qquad k=1,2,\ldots ; \]
where \( T_0 (h) = h^2 \,y(0)^2 \) and
\[ T_{k} = T_{k-1} (h) + \frac{h}{2} \left( y_k + y_{k-1} \right) ,\qquad k=1,2,\ldots ; \]
Finally, we ask matlab to perform all calculations:
h = .001; % step
a = 0;    % initial point
b = 2;    % terminal point
y0 = 50;  % initial value y(a)
n = floor((b-a)/h); % the number of steps (rounded);
y = zeros(n+1,1);   % memory allocation
T = zeros(n,1);     % memory allocation
y(1)=y0;
T(1) = h^2 * y0^2;
y(2) = y0 + 2.3*h*y0 - 0.01*h*y0^2 - 0.1*T(1);
for k= 2:n
  T(k) = T(k - 1) + h*(y(k) + y(k - 1))/2;
  y(k + 1) = y(k) + h*y(k)*(2.3 - 0.01*y(k) - 0.1*T(k));
end
plot(h*(0:n),y(1:n+1))

There are several ways to implement the Euler numerical method for solving initial value problems. We demonstrate its implementation in a series of codes. To start, define the initial point and then the slope function:

Clear[x0,y0,x,y,f]
{x0, y0} = {0, 1}
f[x_, y_] = x^2 + y^2          (* slope function f(x,y) = x^2 + y^2 was chosen for concreteness *)

Next, define the step size:

h = 0.1
Now we define the Euler method itself:
euler[{x_, y_}] = {x + h, y + h*f[x, y]}
Create the table of approximations using Euler's rule:
eilist = NestList[euler, {x0, y0}, 10]

Plot with some options:

plp = ListPlot[eilist]

or
ListPlot[eilist, Joined -> True]
or
ListPlot[eilist, Joined -> True, Mesh -> All]
or
ListPlot[eilist, Filling -> Axis]

Another way is to make a loop.

Clear[y]
y[0]=1; f[x_,y_]=x^2 +y^2 ;
Do[y[n + 1] = y[n] + h f[1+.01 n, y[n]], {n, 0, 99}]
y[10]

Explanation. First of all, it is always important to clear all previous assignment to all the variables that we are going to use, so we have to type: Clear[y]

The basic structure of the loop is:
Do[some expression with n, {n, starting number, end number}]
or with option increment, denoted by k:
Do[some expression with n, {n, starting number, end number, k}]

The function of this Do loop is to repeat the expression, with n taking values from “starting number” to “ending number,” and therefore repeat the expression (1+ (ending number) - (starting number))  many times. For our example, we want to iterate 99 steps, so n will go from 0, 1, 2, ..., until 99, which is 100 steps in total. There is just one technical issue: we must have n as an integer. Therefore we let y[n] denote the actual value of y(x) when x = n*h. This way everything is an integer, and we have our nice do loop.
Finally put y[10] , which is actualy y(0.1), then shift+return, and we have our nice answer.

First, we start with output of our program, which is, perhaps, the most important part of our program design. You'll notice that in the program description, we are told that the "solution that it produces will be returned to the user in the form of a list of points. The user can then do whatever one likes with this output, such as create a graph, or utilize the point estimates for other purposes.

Next, we must decide upon the information that must be passed to the program for it to be able to achieve this. In programming jargon, these values are often referred to as parameters. So what parameters would an Euler program need to know in order to start solving the problem numerically? Quite a few should come to mind:

The slope function f(x, y)
The initial x-value, x0
The initial y-value, y0
The final x-value,
The number of subdivisions to be used when chopping the solution interval into individual jumps or the step size.

To code these parameters into the program, we need to decide on actual variable names for each of them. Let's choose variable names as follows:

f, for the slope function f(x, y)
x0,      for the initial x-value, x0
y0,      for the initial y-value, y0
xn,      for the final x-value
Steps, for the number of subdivisions
h,        for step size

There are many ways to achieve this goal. However, it is natural to have the new subroutine that looks similar to build-in Mathematica's commands. So we want our program might look something like this:

euler[f,{x,x0,xn},{y,y0},Steps]

In order to use Euler's method to generate a numerical solution to an initial value problem of the form:

\[ y' = f(x, y) , \qquad y(x_0 ) = y_0 . \]

We have to decide upon what interval, starting at the initial point x0, we desire to find the solution. We chop this interval into small subdivisions of length h, called step size. Then, using the initial condition \( (x_0,y_0) \) as our starting point, we generate the rest of the solution by using the iterative formulas:

\[ \begin{split} x_{n+1} = x_n + h , \\ y_{n+1} = y_n + h f(x_n, y_n) \end{split} \]
to find the coordinates of the points in our numerical solution. The algorithm is terminated when the right end of the desired interval has been reached.

If you would like to use a built-in Euler’s method, unluckily there is none. However, we can define it ourself. Simply copy the following command line by line:

euler[f_, {x_, x0_, xn_}, {y_, y0_}, Steps_] :=
Block[{ xold = x0, yold = y0, sollist = {{x0, y0}}, h },
h = N[(xn - x0) /Steps];            (* or n = (xn - x0) /Steps//N; *)
Do[ xnew = xold + h;
ynew = yold + h * (f /. {x -> xold, y -> yold});
sollist = Append[sollist, {xnew, ynew}];
xold = xnew;
yold = ynew,
{Steps}
];
Return[sollist]
]

Now we have our euler function:   euler[f(x,y), {x,x0,x1},{y,y0},steps]
Then this script will solve the differential equation y’=f(x,y), subject to the initial condition y(x0)=y0, and generate all values between x0 and x1. The number of steps for the Euler’s method is specified with steps.
To solve the same problem as above, we simply need to input:

euler[1/(3*x - 2*y + 1), {x, 0, 0.4}, {y, 1}, 4]
Out[2]= {{0, 1}, {0.1, 0.9}, {0.2, 0.7}, {0.3, 1.2}, {0.4, 1.}}

Another example:

euler[y+x, {x, 0, 1}, {y, 1}, 10]
Out[3]= {{0, 1}, {0.1, 1.1}, {0.2, 1.22}, {0.3, 1.362}, {0.4, 1.5282}, {0.5,
1.72102}, {0.6, 1.94312}, {0.7, 2.19743}, {0.8, 2.48718}, {0.9, 2.8159}, {1., 3.18748}}

As you can see, the output is really a big table of values, and you can really just read the last one off to get y(1). If the only one final value is needed, then Return command should be replaced with
Return[sollist[[Steps + 1]]]]

We can plot the output as follows

ListPlot[euler[y + x, {x, 0, 1}, {y, 1}, 10]]

Or we can plot it as

aa = euler[y + x, {x, 0, 1}, {y, 1}, 10]
Out[4]= {{0, 1}, {0.1, 1.1}, {0.2, 1.22}, {0.3, 1.362}, {0.4, 1.5282}, {0.5,
1.72102}, {0.6, 1.94312}, {0.7, 2.19743}, {0.8, 2.48718}, {0.9,
2.8159}, {1., 3.18748}}
ListPlot[aa, AxesLabel -> {"x", "y"}, PlotStyle -> {PointSize[0.015]}]


The following code, which uses a slightly different programming paradigm, implements Euler's method for a system of differential equations:

euler[F_, a_, Y0_, b_, Steps_] :=
Module[{t, Y, h = (b - a)/Steps//N, i},
t[0] = a; Y[0] = Y0;
Do[
t[i] = a + i h;
Y[i] = Y[i - 1] + h F[t[i - 1], Y[i - 1]],
{i, 1, n}
];
Table[{t[i], Y[i]}, {i, 0, n}]
]

And the usage message is:

euler::usage = "euler[F, t0, Y0, b, n] gives the numerical solution to {Y' == F[t, Y], Y[t0] == Y0} over the interval\n
[t0, b] by the n-step Euler's method. The result is in the form of a table of {t, Y} pairs."

Note that this function uses an exact increment h rather than converting it explicitly to numeric form using N. Of course you can readily hange that or accomplish the corresponding thing by giving an approximate number for one of the endpoints.

Next we plot the points by joining them with a curve:

a = ListPlot[euler[f, 0, 1, 3, 30], Joined -> True]

Another way without writing a subroutine:

f[x_, y_] := y^2 - 3*x^2;
x0 = 0;
y0 = 1;
xend = 1.1;
steps = 5;
h = (xend - x0)/steps // N;
x = x0;
y = y0;
eulerlist = {{x, y}};
For[i = 1, i <= steps, y = f[x, y]*h + y;
x = x + h;
eulerlist = Append[eulerlist, {x, y}]; i++]
Print[eulerlist]

The results can also be visualized by connecting the points:

a = ListPlot[eulerlist, Joined -> True, Epilog -> {PointSize[0.02], Map[Point, eulerlist]}];
s = NDSolve[{u'[t] == f[x, u[t]], u[0] == 1}, u[t], {t, 0, 1.1}];
b = Plot[Evaluate[u[x] /. s], {x, 0, 1.1}, PlotStyle -> RGBColor[1, 0, 0]];
Show[a,b]

Next, we demonstrate application of Function command:

EulerODE[f_ /; Head[f] == Function, {t0_, y0_}, t1_, n_] :=
Module[{h = (t1 - t0)/n // N, tt, yy},
tt[0] = t0; yy[0] = y0;
tt[k_ /; 0 < k < n] := tt[k] = tt[k - 1] + h;
yy[k_ /; 0 < k < n] :=
yy[k] = yy[k - 1] + h f[tt[k - 1], yy[k - 1]];
Table[{tt[k], yy[k]}, {k, 0, n - 1}]
];

ty = EulerODE[Function[{t, y}, y^2/t/2 + y^2/t^1.5], {1, 1}, 2, 100];
Plot[Interpolation[ty][t], {t, 1, 2}]

Here is another steamline of the Euler method (based on Mathematica function).

Clear[ x,y,h,i]
f[x_, y_] := y^2 - 3*x^2
x[1] = 0;     y[1] = 1;
h = 0.1;
Do[ { x[i+1] = x[i]+h,
          y[i+1] = y[i] +h*f[x[i],y[i]] } , {i,1,20} ]

Now we can plot our results, making sure we only refer to values of x and y that we have defined:

Clear[pairs]
pairs = Table[{x[i], y[i]}, {i, 1, 11}]
ListPlot[pairs, Joined -> True]
Out[23]= {{0, 1}, {0.1, 1.1}, {0.2, 1.218}, {0.3, 1.35435}, {0.4, 1.51078}, {0.5, 1.69102}, {0.6, 1.90198}, {0.7, 2.15573}, {0.8, 2.47345}, {0.9, 2.89325}, {1., 3.48734}}
Finally, we plot Euler approximations along with the actual solution:
soln = DSolve[{f[x, y[x]] == y'[x], y[0] == -1}, y, x]
plot1 = Plot[Evaluate[y[x] /. soln], {x, 0, 2.5}, PlotStyle -> {{Thickness[0.008], RGBColor[1, 0, 0]}}];
Show[plot1, plot5]

Now we are going to repeat the problem, but using Mathematica lists format instead. The idea is still the same---we set a new x value, compute the corresponding y value, and save the two quantities in two lists. In particular, in a list we are simply storing values; In previously presented code, we were actually storing rules, which are more complicated to store and evaluate. We start out with just one value in each list (the initial conditions). Then we use the Append command to add our next pair of values to the list.

Clear[ x,y,h,i]
x = {0.0 };
y = {2.0 };
h = 0.1;
Do [ { x = Append[x,x[[i]]+h],
           y = Append[y,y[[i]]+h*f[x[[i]],y[[i]] ] ],
         } , {i,1,20} ]

Now we can plot our results:

Clear[ pairs ]
pairs = Table[ {x[[i]], y[[i]] }, {i,1,21} ]
ListPlot [ pairs, Joined -> True ]

Another option:
f[x_, y_] := Exp[2*x - y]
h = 0.1
K = IntegerPart[1.2/h]  (* if the value at x=1.2 is needed *)
y[0] = 1                      (* initial condition *)


Euler method:
Do[y[n + 1] = y[n] + h*f[n*h, y[n]], {n, 0, K}]
y[K]

Euler's method is first-order accurate, so that errors occur one order higher starting at powers of \( h^2 . \)
Euler's method is implemented in NDSolve as ExplicitEuler:

NDSolve[{y'[t] == t^2 - y[t], y[0] == 1}, y[t], {t, 0, 2},
Method -> "ExplicitEuler", "StartingStepSize" -> 1/10]

 

  1. Atkinson, K., Han, W., Stewart, D., Numerical Solutions of Ordinary Differential Equations, Wiley, 2009.
  2. Atkinson, K., Han, W., Stewart, D., ODE Programs
  3. Hataue, I., Analysis of ghost numerical solutions of differential equation caused by nonlinear instability, 2012 https://doi.org/10.2514/3.12350
  4. Numerical methods for Differential Equations.
  5. Yamaguti, M., Ushiki, S., Chaos in numerical analysis of ordinary differential equations, Physica D: Nonlinear Phenomena Volume 3, Issue 3, August 1981, Pages 618-626; https://doi.org/10.1016/0167-2789(81)90044-0