Clear[a,b,x,y,f,h,xx,yy,"k*"]
f[x_,y_]:=Sqrt[x+y]+ y*Cos[x y];
a=1;b=2;h=0.005;
n=((b-a)/h);
xx=Table[a+i*h,{i,0,n}];
yy=Table[0,{i,0,n}]; (* декларираме структурата на yy: мрежовите стойности на функцията, които търсим *)
yy[[1]]=1;
For[i=1,i ≤ n,i++,
(*body*)
k1=h*f[xx[[i]],yy[[i]]];
k2=h*f[xx[[i]]+(1/2)*h,yy[[i]]+(1/2)*k1];
k3=h*f[xx[[i]]+(1/2)*h,yy[[i]]+(1/2)*k2];
k4=h*f[xx[[i]]+h,yy[[i]]+k3];
yy[[i+1]]=yy[[i]]+(1/6)*(k1+2 k2 +2 k3 +k4)
(* end For *)]
TableForm[SetPrecision[{xx,yy,f[xx,yy]},10],TableHeadings->{{"x","y","f(x,y)"},Table[i,{i,0,n}]}]
points=Table[{xx[[i]], yy[[i]]},{i,1,n+1}];
ListPlot[points]

точност на резултата:

(h^4)