Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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]-Cx[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]); 

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 suppose to fluctuate around the equilibrium point but in the code I have so far, this doesn't happen.

share|improve this question
Welcome. Why don't you show what you have done so far? – user21 Apr 18 '13 at 4:48
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 – andre Apr 18 '13 at 7:53
up vote 14 down vote accepted

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 this 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.

share|improve this answer
I appreciate your help, I have achieved my simulation. With the huge help of your answer! – Isaias Cordova 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

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

yl = MapThread[Append, {yf["Grid"], yf["ValuesOnGrid"]}]
   {{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]]

plot of approximate solution from RK

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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