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[{},
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);

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;

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?

  • 2
    $\begingroup$ there are many problems in your code. First you do not find k4. Second, you need ; between end of the For loops and the Return[] statement. Also you do not need an explicit Return. Also the call is wrong. Why not just RK[0, 2, 2, 5, f[x, x]] ? Try to correct these first and see. btw, there is lots of RK4 code in this forum, many questions were asked about it before. If you google, you'll find examples. $\endgroup$
    – Nasser
    Commented Apr 22, 2020 at 17:26
  • $\begingroup$ The k4 I just messed up in post...its correct in my code. Ill fix the others, thank you. $\endgroup$ Commented Apr 22, 2020 at 17:28
  • $\begingroup$ I am not sure how f[x,x] is equivalent...or correct. Wont that just solve x - x^4? $\endgroup$ Commented Apr 22, 2020 at 17:30
  • 1
    $\begingroup$ You also failed to localize any variables within Module. $\endgroup$
    – Szabolcs
    Commented Apr 22, 2020 at 17:34

2 Answers 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

Mathematica graphics

  • $\begingroup$ Nice fix on the output... the first one destroyed my browser! Thanks for being so helpful! $\endgroup$ Commented Apr 22, 2020 at 17:47

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))$.

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 = 
  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}},

 Transpose[{res[[All, 1]], res[[All, 2, 1]]}]
  • $\begingroup$ Thanks for the spunked up version. I am a Mathematica Ignoramus, so I am gonna stick with Nasser's fix on my code. Yours looks cool, but its beyond me! Thanks for being so helpful! $\endgroup$ Commented Apr 22, 2020 at 17:43
  • $\begingroup$ @ragnvaldr.js Recommended reading: mathematica.stackexchange.com/q/134609/12 $\endgroup$
    – Szabolcs
    Commented Apr 22, 2020 at 17:44
  • $\begingroup$ Thanks... I had no clue that the For loop should be avoided in Mathematica. For what I am doing, this will work, but Ill try and work through your solution for future use... Thanks! $\endgroup$ Commented Apr 22, 2020 at 17:49
  • 1
    $\begingroup$ I'm sure it is... but I'm an amatuer Java Developer so the For loop is just second nature to me. I guess I am confused as the res = ... stuff. The added info helps... The syntax is still just unfamiliar. $\endgroup$ Commented Apr 22, 2020 at 17:56
  • 1
    $\begingroup$ @Szabolcs Your code is perfect (+1), see my answer on mathematica.stackexchange.com/questions/221848/… $\endgroup$ Commented May 14, 2020 at 22:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.