Numerical Methods of Solving 1st order ODEs:  Euler's Method

The initial value problem and the exact solution

Consider the DE y ' = f(x, y) = (e^(2x)) + y with the IC y[0]=5. 
(If our initial condition were y(5)=10, we'd have x[0]=5, y[0]=10.)
This is easy to find the exact solution, so we do not need to use numerical techniques.  
We use it in order to illustrate how well (or poorly) the method works.  This method is simple but not necessarily accurate.

Here's the functional part before the initial condition.  It is this family of functions that we can visualize in the direction field.

Clear[x,y]
DSolve[y'[x]==(E^(2x))+ y[x],y[x],x]

The unique solution to the initial value problem:

DSolve[{y'[x]==(E^(2x))+ y[x],y[0]==5},y[x],x]

soln[x_]:=(E^x) (4+(E^x))

Clear[p1]
p1=Plot[(E^x) (4+(E^x)),{x,0,0.6}]

six manual iterations with step size h=.1

The initial set-up

Clear[n,x,y,h,slope, f]

h=.1;

n=0;
x[n]=0;
y[n]=5;

f[x_,y_]:=(E^(2x))+ y

The iterations

We begin with steps repeated (copied and pasted) to ease us into using a FOR loop.

n=n+1
x[n]=x[n-1] + h
slope  = f[ x[n-1], y[n-1] ]
y[n] = y[n-1] + h * slope

n=n+1
x[n]=x[n-1] + h
slope  = f[ x[n-1], y[n-1] ]
y[n] = y[n-1] + h * slope

n=n+1
x[n]=x[n-1] + h
slope  = f[ x[n-1], y[n-1] ]
y[n] = y[n-1] + h * slope

n=n+1
x[n]=x[n-1] + h
slope  = f[ x[n-1], y[n-1] ]
y[n] = y[n-1] + h * slope

n=n+1
x[n]=x[n-1] + h
slope  = f[ x[n-1], y[n-1] ]
y[n] = y[n-1] + h * slope

n=n+1
x[n]=x[n-1] + h
slope  = f[ x[n-1], y[n-1] ]
y[n] = y[n-1] + h * slope

The results in table form

n now represents the number of iterations we made, so we use another variable, i, to march through the values from 0 to n.

TableForm[Table[{i,x[i],y[i]},{i,0,n}]]

Compare with the known solution

Relative Error

Because this is an instructional example, we calculate the actual value and relative error.

TableForm[Table[{y[i],soln[x[i]], Abs[(soln[x[i]] - y[i])/y[i]]},{i,0,n}]]

Clear[pointlist]
pointlist=Table[{x[i],y[i]},{i,0,n}]

Clear[p2]
p2=ListPlot[pointlist,Joined->True,PlotStyle->{Dashing[{0.15,0.05}],Hue[0.05]}]

As you can see, our numerical solution is very close for just a little while.  In just 6 steps, though, it's pretty far off.

Show[{p1,p2}]

Twelve iterations with h=.05

As you can see, our numerical solution is much better with a smaller step size.

The initial set-up

Clear[n,x,y,h,slope]
h=.05;

x[0]=0;
y[0]=5;

The iterations in a loop

For[n=0,n≤12,
n=n+1;
x[n]=x[n-1] + h;
slope  = f[ x[n-1], y[n-1] ];
y[n] = y[n-1] + h * slope]

The results in table form

TableForm[Table[{i,x[i],y[i]},{i,0,n}]]

Compare with the known solution

Clear[pointlist]
pointlist=Table[{x[i],y[i]},{i,0,n}]

TableForm[Table[{y[i],soln[x[i]], Abs[(soln[x[i]] - y[i])/y[i]]},{i,0,n}]]

Clear[p3]
p3=ListPlot[pointlist]

Show[{p1,p3}]

Show[{p1,p2,p3}]

Compare with direction field

Even when the DE is hard or impossible to solve, we can compare it to the direction field.

p4=VectorPlot[{1,f[x,y]},{x,0.0,2.0},{y,4.0,11.0},Axes->True]

Show[{p1,p2,p3,p4},Axes->True]