I am trying to solve differential equations numerically, so I am trying to write a 4th -order Runga-Kutta program for Mathematica (I know NDSolve does this, but I want to do my own). I ran into some trouble though, as my program just loops infinitely.

RK[a_,b_,y0_,n_,f_]:= Module[{},
h=(b-a)/n;
X = Table[a+k*h, {k,0,n} ];
Y = Table[y0, {k,0,n} ];
For[j=1, j<n, j++,
k1 = f[X[[j]],Y[[j]]];
k2 = f[X[[j]]+(h/2),Y[[j]]+h*(k1/2)];
k3 = f[X[[j]]+(h/2),Y[[j]]+h*(k2/2)];
k4 = f[X[[j+1]],Y[[j]]+h*k3];
Y[[j+1]]= Y[[j]]+(h/6)(k1+2*k2+2*k3+k4);
];
Return[Transpose[{X,Y}]];
];

I don't think my issue is with the algortithm though... I think its with my definition of the differential equation. I was honestly pretty lost on how I do this, but this is what I came up with:

f[x_,y_] = y - (x^2)(y)^2;
RK[0,10,2,50,f[x,Function[x,y[x]]]]

I tried defining it as a function of two variables... but I think I might have done some thing wrong.

If this is wrong...how do I define a differential equation as a function of two variables?