Numerical Methods of Solving 1st order ODEs:  Runga-Kutta Method

The initial value problem

Consider the DE y ' = f(x, y) =  x - y, y(0)=1. 
This is easy to find the exact solution, so we do not need to use numerical techniques.

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] == x-y[x], y[x], x]

The unique solution to the initial value problem:

DSolve[{y'[x] == x-y[x], y[0]==1}, y[x], x]

soln[x_]:=(E^(-x)) (2-(E^x)+(E^x) x)
Clear[p1]
p1=Plot[soln[x],{x,0.0,5.0}]

iterations with  step size h=.05

Clear[n,h,x,y,f,m1,m2,m3,m4,numIterations]

x[0]=0.0;
y[0]=1.0;

h=.05;

f[x_,y_]:=x - y

numIterations = 100;

For[n=0,n≤numIterations,

m1 = f[ x[n], y[n] ];
m2  = f[ x[n] + h/2 , y[n] + m1*h/2];
m3  = f[ x[n] + h/2 , y[n] + m2*h/2];
m4  = f[ x[n] + h , y[n] + m3*h];

n=n+1;
x[n]=x[n-1] + h;
y[n] = y[n-1] + (m1+2*m2+2*m3+m4)*h/6;     ]
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]}]

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])/soln[x[i]] ] },{i,1,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]}];
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.

Clear[p3]
p3=VectorPlot[{1,f[x,y]},{x,-0.5,5.5},{y,-0.5,4.5}];
Show[{p1,p2,p3},Axes->True]

Show[{p1,p2,p3},Axes->True,PlotRange->{{-0.5,5.5},{-0.5,4.5}},AxesOrigin->{0,0}]