Numerical Methods of Solving 1st order ODEs:  Improved Euler's Method or Heun 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 not as simple as the first Euler method but is more 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 iterations with step size h=.1

The initial set-up

Clear[n,x,y,h,f,slope1,slope2]

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

h=.1;

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

For[n=0,n≤6,
n=n+1;
x[n]=x[n-1] + h;
slope1  = f[ x[n-1], y[n-1] ];
slope2  = f[ x[n], y[n-1] + h * slope1 ];
y[n] = y[n-1] + h * (slope1 + slope2)/2]

The results in table form

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

Compare with the known solution

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

The numerical solution is excellent on a larger, though still small, interval.

Show[{p1,p2}]

Compare with direction field

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

Needs["VectorFieldPlots`"];
p3=VectorFieldPlots`VectorFieldPlot[{1,f[x,y]},{x,0.0,2.0},{y,4.0,11.0}]

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

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

solve on a larger interval

Clear[n,x,y,h,slope1,slope2, numIterations,p3, p4,pointlist]

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

h=.1;

numIterations = 50;

For[n=0,n≤numIterations,
n=n+1;
x[n]=x[n-1] + h;
slope1  = f[ x[n-1], y[n-1] ];
slope2  = f[ x[n], y[n-1] + h * slope1 ];
y[n] = y[n-1] + h * (slope1 + slope2)/2]
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}]]

p3=Plot[(E^x) (4+(E^x)),{x,0,numIterations * h}];
p4=ListPlot[pointlist,Joined->True,PlotStyle->{Dashing[{0.15,0.05}],Hue[0.05]}];
Show[{p3,p4}]