Mathematica Runge-Kutta

runge_kutta.nb

(* One step Runge-Kutta                                                           *)
(*  fn    : Definition of the function                                            *)
(*  init  : List of initial values, dimensions have to agree with the input of fn *)
(*  step  : Step size for the numerical simulation                                *)
(* Usage:                                                                         *)
(*  One step:                                                                     *)
(*    RungeKutta[{#[[2]], #[[1]]^2 - 2#[[2]]}&,{.1, -.5}, .1]                     *)
(*  Multiple steps:                                                               *)
(*    NestList[RungeKutta[{#[[2]], #[[1]]^2 - 2#[[2]]}&, #, 0.1]& , {.1, -.5}, 5] *)

RungeKutta[fn_, init_, step_] := Module[
  {k1, k2, k3, k4},
  k1 = step fn[init];
  k2 = step fn[init + k1 / 2];
  k3 = step fn[init + k2 / 2];
  k4 = step fn[init + k3];
  init + Total[{k1, 2 k2, 2 k3, k4}]/6
]

van_der_pol_example.nb

stepSize = 0.0004;
initialValues = {2.12, -3.1};
vanDerPolParam = 5;
nrIterations = 100000;

vanDerPol := {#[[2]], #2 (1 - #[[1]]^2) #[[2]] - #[[1]]} &

RungeKutta[fn_, init_, step_] := Module[
  {k1, k2, k3, k4},
  k1 = step fn[init];
  k2 = step fn[init + k1 / 2];
  k3 = step fn[init + k2 / 2];
  k4 = step fn[init + k3];
  init + Total[{k1, 2 k2, 2 k3, k4}]/6
]

ListPlot[NestList[
  RungeKutta[(vanDerPol[#, vanDerPolParam]) &, #, stepSize] & , 
  initialValues, nrIterations],
 Frame -> True,
 Axes -> None,
 AspectRatio -> 1,
 PlotRange -> All]1