# Solving a system of ODEs with the Runge-Kutta method

I'm trying to solve a system of ODEs using a fourth-order Runge-Kutta method. I have to recreate certain results to obtain my degree. But I'm a beginner at Mathematica programming and with the Runge-Kutta method as well.

{A = 0.30, B = 1, C = 40, D = 1, E = 0.75, F = 0.11,
r = 2.5, a = 2, e = 0.475, g = 2, d = 0.03, n = 0.01, p = -0.00005}

x'[t]/x[t] = (A + B x[t] - C x[t]^2 - F)/(D + E^(r y[t]));
y'[t]/y[t] = g (((a s[t] x[t] k[t])/m[t]) - e) (1 - y[t]);
m'[t]/m[t] = n;
k'[t]/k[t] = x[t] - d;
s'[t]/s[t] = -p;


I'd appreciate any kind of help. For over a month now, I've tried to solve this system myself but have only gotten bad results. This model is supposed to fluctuate around the equilibrium point, but in the code I have so far, this doesn't happen.

• Please don't use single capital letters as variables: both C and E are protected system symbols, and no value should be assigned to them. Also note that multiplication requires at least a space, thus Cx is a symbol while C x is C times x Commented Apr 18, 2013 at 7:50
• Runge-Kutta 4 is described in the documentation as a example here Commented Apr 18, 2013 at 7:53

Here is a functional approach. The following will give you one step of the Runge-Kutta formula:

RungeKutta[func_List, yinit_List, y_List, step_] :=
Module[{k1, k2, k3, k4},
k1 = step N[func /. MapThread[Rule, {y, yinit}]];
k2 = step N[func /. MapThread[Rule, {y, k1/2 + yinit}]];
k3 = step N[func /. MapThread[Rule, {y, k2/2 + yinit}]];
k4 = step N[func /. MapThread[Rule, {y, k3 + yinit}]];
yinit + Total[{k1, 2 k2, 2 k3, k4}]/6]


Here, func is a list of functions, yinit a list of initial values, y a list of function variables (in your case that will be {x, y, m, k, s}, and step is the step size of the numerical simulation. You can then use something like NestList to iterate it many times like this:

NestList[RungeKutta[func, #, y, step] &, N[yinit], Round[t/step]]


Here, t is the maximum value of your independent variable. You can also include conditionals to check that the lists are of equal length.

• I appreciate your help, I have achieved my simulation. With the huge help of your answer!
– user6983
Commented Apr 19, 2013 at 0:34
• @IsaiasCordova, glad I could help. Commented Apr 19, 2013 at 0:51
• Er… I think it's better to mention that the y for RungeKutta[] usually include t for more general cases. (And the counterpart of t is 1 in yinit.) Commented Jun 26, 2013 at 10:51
– TMH
Commented Apr 9, 2014 at 14:11
• @xzczd, yes, that version as written in this answer is intended for autonomous equations; Maeder's book implemented it this way initially, but then also gave a version for non-autonomous DEs. Commented Mar 28, 2018 at 0:49

andre has linked you to the method plug-in framework in the comments, but there is a more direct way to implement classical Runge-Kutta, just by supplying its Butcher tableau to "ExplicitRungeKutta". Here's an adaptation from the docs:

ClassicalRungeKuttaCoefficients[4, prec_] := With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}},
bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, N[{amat, bvec, cvec}, prec]]

{xf, yf} = {x, y} /. First @
NDSolve[{x'[t] == -y[t], y'[t] == x[t], x[0] == 1, y[0] == 0},
{x, y}, {t, 0, 6},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
"Coefficients" -> ClassicalRungeKuttaCoefficients},
StartingStepSize -> 1/2];


To obtain the values computed by RK, you can then do this:

xl = MapThread[Append, {xf["Grid"], xf["ValuesOnGrid"]}]
{{0., 1.}, {0.5, 0.877604}, {1., 0.540588}, {1.5, 0.0714256}, {2., -0.415108},
{2.5, -0.800012}, {3., -0.989166}, {3.5, -0.936349}, {4., -0.65453}, {4.5, -0.212684},
{5., 0.281088}, {5.5, 0.706007}, {6., 0.95816}}

{{0., 0.}, {0.5, 0.479167}, {1., 0.841037}, {1.5, 0.99713}, {2., 0.90931},
{2.5, 0.599108}, {3., 0.142441}, {3.5, -0.348969}, {4., -0.754924}, {4.5, -0.976153},
{5., -0.958587}, {5.5, -0.706572}, {6., -0.281796}}


Plot them:

ListLinePlot[{xl, yl}, Mesh -> All, MeshStyle -> PointSize[Medium]]


• Hey, I have a question. why do they call the vector on the inner part of the butchear tableau amac and not avec? I don't quite see why they call it that. Not a big issue, but just conufsed by it. It is because it's a matrix and not a vector? and if so, why doesn't it have for example, as the entries, {{1/2,0,0},{0,1/2,0},{0,0,1}}? Is it because we add the factors in as necessary, so exclude and 0s after the last entry for that step of the algorithm? Commented Oct 13, 2019 at 18:00
• @Shinaolord, in fact, the coefficient matrix is a strictly lower triangular matrix, so the actual complete matrix is something like {{0, 0, 0, 0}, {1/2, 0, 0, 0}, {0, 1/2, 0, 0}, {0, 0, 1, 0}}. For explicit RK methods, by convention, the upper triangular entries are omitted; for implicit RK, methods, however, the full matrix will usually be shown in the Butcher tableau. Commented Nov 13, 2019 at 20:19

First of all, I'd like to emphasize that, if you just want to solve an initial value problem (IVP) of ordinary differential equation (ODE) or ODE system, please use NDSolve. Forget about the classical Runge-Kutta (classical RK, RK4), it's mainly for pedagogical purposes nowadays, and not comparable to the default method of NDSolve at all. If default method of NDSolve fails to solve your IVP, RK4 won't save your day, either. Actually I've never seen a knotty problem solved by switching from default NDSolve method to RK4 since I began using Mathematica 10 years ago.

If you need RK4 for certain special reason (e.g. to exactly reproduce numeric result obtained with RK4 in literature, to ensure the result can be easily reproduced by others, etc.), setting proper options for NDSolve is the way to go. Please refer to J.M.'s answer for more info. Here I'll provide a general-purposed package for setting up explicit Runge-Kutta method with no adaptive step control, which is implemented without the help of NDSolve, just for fun:

BeginPackage["primaryRK"];
rkstepgenerator; rhsfuncgenerator; rk4stepfunc;
Begin["private"];

rule = {semi -> CompoundExpression, set -> Set,
module -> Module, compile -> Compile};
SetAttributes[semi, Flat];
f[{}, _] := {};

rkstepgenerator[a_, b_, c_] := (s = Length@b;
k = ToExpression["k" <> ToString@#] & /@ Range@s;
Function @@ {f,
compile[{t, {y, _Real, 1}, dt},
module[k, semi[
set[k, Prepend[Thread@f[t + c dt, y + dt k . PadRight[#, s] & /@ a],
f[t, y]]],
y + dt b . k]]]} /. rule)

a = {{1/2}, {0, 1/2}, {0, 0, 1}};
b = {1/6, 1/3, 1/3, 1/6};
c = {1/2, 1/2, 1};

rk4stepfunc = rkstepgenerator[a, b, c];

rhsfuncgenerator[t_, ylst_][expr_] :=
Function @@ {{t, yvec}, module[ylst, semi[ylst~set~yvec, expr]]} /. rule

End[];
EndPackage[]


Since RK4 is so frequently requested, a function rk4stepfunc is included in the package. Let me use the IVP

\begin{aligned}x'(t)&=-y(t)\\ y'(t)&=x(t)+\sin(3t)\\ x(0)&=1\\ y(0)&=0\end{aligned}

to illustrate its possible usage:

heads = {x, y};
eqs = {x'[t] == -y[t], y'[t] == x[t] + Sin[3 t]};
iclst = {1, 0};
domain = {t0, tend} = {0, 6};

(* Analytic solution for comparison: *)
ref = Plot[
vars /. DSolve[{eqs, vars == iclst /. t -> t0}, vars, t] //
Evaluate, {t, t0, tend}];

expr =
Solve[eqs, D[vars, t]][[1, All, -1]] /.
(h : Alternatives @@ heads)[t] :> h
(* {-y, x + Sin[3 t]} *)

Internal\$ContextMarks = False; (* Default value is Automatic *)

(* Function[{t, yvec}, Module[{x, y}, {x, y} = yvec; {-y, x + Sin[3 t]}]] *)

step = rk4stepfunc@rhsfunc;

dt = 1/30;
tlst = Range[t0, tend, dt];

sollst =
FoldList[step[#2[[1]], #, #2[[2]]] &, iclst,
Transpose@{Most@tlst, Differences@tlst}]; // AbsoluteTiming
(* {0.0013601, Null} *)



You can of course use the package to set up other explicit Runge-Kutta method with no error estimation e.g. explicit Euler. Just for fun, let me show another possible usage of the package:

eulerstepfunc = rkstepgenerator[{{}}, {1}, {}]
(* Function[f,
Compile[{t, {y, _Real, 1}, dt},
Module[{k1}, k1 = f[t, y]; dt k1 + y]]] *)

step = eulerstepfunc@rhsfunc;
dt = 1/50;
ti = t0 - dt;
sollst = NestList[
step[ti = ti + dt, #, dt] &, iclst, #2 - # & @@ domain/dt //
Round]; // AbsoluteTiming
(* {0.0007594, Null} *)
ListPlot[sollst\[Transpose], DataRange -> domain, PlotStyle -> Red]~Show~ref


Heun's method (Once again, I'll show a different usage of the package for fun):

heunstepfunc =
With[{α = 1},
rkstepgenerator[{{α}}, {1 - 1/(2 α), 1/(2 α)}, {α}]]
(* Function[f,
Compile[{t, {y, _Real, 1}, dt},
Module[{k1, k2}, k1 = f[t, y];
k2 = f[dt + t, dt k1 + y];
dt (k1/2 + k2/2) + y]]] *)

step = heunstepfunc@rhsfunc;
sollst = With[{dt = 1/20, t0 = N@t0, step = step, iclst = N@iclst},
With[{n = 1 - Subtract @@ domain/dt},
Compile[{}, Module[{lst = Table[iclst, {n}], ti = t0},
Do[lst[[i]] = step[ti, lst[[i - 1]], dt];
ti = ti + dt, {i, 2, n}];
lst]]]][]; // AbsoluteTiming

Show[sollst\[Transpose] //
ListPlot[#, DataRange -> domain, ColorFunction -> "Rainbow"] &, ref]


• Very nice indeed! :) Commented Oct 1, 2022 at 13:59
• Thanks for emphasizing NDSolve. There are so many questions here about RK4. It seems many supervisors have zero idea about the capabilities of modern computer algebra systems. It's like teaching driving cars with manual transmission when everyone is already driving Tesla. Commented Oct 1, 2022 at 14:13
• @yarchik "It seems many supervisors have zero idea..." — I've seen the Euler method used fairly often. Like using a hand crank to start your car. :) Commented Oct 1, 2022 at 15:46