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

The initial value problem

Consider the DE y ' = f(x, y) = sin(x) + cos(y) with the IC y[0]=0. 
This is hard to find the exact solution, so we 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] == Sin[x] + Cos[y[x]], y[x], x]

Solve::"ifun": "Inverse functions are being used by \!\(Solve\), so some solutions may not be found; use Reduce for complete solution information. \!\(\*ButtonBox[\"»\", ButtonStyle->\"Link\", ButtonFrame->None, ButtonData:>\"paclet:ref/message/Solve/ifun\", ButtonNote -> \"Solve::ifun\"]\)"

$Aborted

iterations with  step size h=.1

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

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

h=.1;

f[x_,y_]:=Sin[x]+Cos[y]

numIterations = 50;

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 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.0,5.5},{y,0.0,3.5}];
Show[{p2,p3},Axes->True]