Runge-Kutta implemented on Mathematica

Question

I am trying to solve differential equations numerically, so I am trying to write a 4th -order Runge-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?

Answer 1

Nasser already pointed out many mistakes, so I won't go into those.

NestList would allow for a much cleaner implementation.

Below, RK4step[f,h] denotes a function which takes a pair of $\{t,y(t)\}$ values, and produces the next one at $t+h$, assuming that $y'(t) = f(t, y(t))$.

ClearAll[RK4step]
RK4step[f_, h_][{t_, y_}] :=
 Module[{k1, k2, k3, k4},
  k1 = f[t,       y];
  k2 = f[t + h/2, y + h k1/2];
  k3 = f[t + h/2, y + h k2/2];
  k4 = f[t + h,   y + h k3];
  {t + h, y + h/6*(k1 + 2 k2 + 2 k3 + k4)}
 ]

We can use NestList to take a starting pair $\{t_0, y(t_0)\}$, and repeatedly propagate the time using RK4step.

res = 
 NestList[
  RK4step[-#2 &, 0.1], (* #2 & is short for f where f[t_, y_] := -y, look up Function *)
  {0.0, 1.0}, (* this is {t0, y(t0)} *)
  100 (* compute this many steps *)
 ]

ListPlot[res, PlotRange -> All]

More complex example, a harmonic oscillator:

f[t_, {x_, v_}] := {v, -x}

res = NestList[
   RK4step[f, 0.1],
   {0.0, {1.0, 0.0}},
   100 
   ];

ListPlot[
 Transpose[{res[[All, 1]], res[[All, 2, 1]]}]
 ]

Answer 2

This works for me

RK[a_, b_, y0_, n_, f_] := Module[{X, Y, j, k1, k2, k3, k4, h},
  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);
   ];

   Transpose[{X, Y}]
  ];

f[x_, y_] := y - (x^2) (y)^2;
RK[0, 2, 2, 5, f] // N