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:

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 𝑎 = x₀ is the starting point) into m equal subintervals and select the mesh/grid points

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 = x_n.
   \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 y_n 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 x_n and x_n+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 y_n, we again obtain the Euler formula.


3. We may assume the possibility of expanding the solution in a Taylor series around the point x_n:
   \[ 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:    y_n+1 = y_n + h f(x_n, y_n).

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]

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 (x_n+1, y_n+1) according to the Euler rule:    y_n+1 = y_n + h f(x_n, y_n). Therefore, every next Euler point y_n+1 is determined according to the tangent line starting from the previous grid point (x_n,y_n) with slope f(x_n,y_n); recall that the equation of the tangent line is y = y_n + f(x_n,y_n) (x - x_n).

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];

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))

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

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

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])

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’]) 

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

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.

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

• 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

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)

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

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

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

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:

Next, define the step size:

Plot with some options:

or

Another way is to make a loop.

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, x₀
The initial y-value, y₀
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)
x₀,      for the initial x-value, x₀
y₀,      for the initial y-value, y₀
x_n,      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:

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

We have to decide upon what interval, starting at the initial point x₀, 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:

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:

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 x₀ 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:

Another example:

We can plot the output as follows

Or we can plot it as

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

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:

Another way without writing a subroutine:

The results can also be visualized by connecting the points:

Next, we demonstrate application of Function command:

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

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

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.

Now we can plot our results:

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:

figure(1)
plot(x, y, 'Color', [1 0 0]) %blue line
hold on
plot(x, z, 'Color', [0 1 0]) %green line