# 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 – István Zachar Apr 18 '13 at 7:50
• Runge-Kutta 4 is described in the documentation as a example here – andre314 Apr 18 '13 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 Apr 19 '13 at 0:34
• @IsaiasCordova, glad I could help. – RunnyKine Apr 19 '13 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.) – xzczd Jun 26 '13 at 10:51
• @IsaiasCordova Can you please post the complete answer here? – TMH Apr 9 '14 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. – J. M.'s torpor Mar 28 '18 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}}


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

• @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. – J. M.'s torpor Nov 13 '19 at 20:19