I am a beginner at Mathematica programming and with the Runge-Kutta method as well. I'm trying to solve a system of coupled ODEs using a 4th-order Runge-Kutta method for my project work.

I have solved it by NDSolve, but I want to solve this by 4th-order Runge-Kutta method. Here is my problem:

Γ = 1.4    
k = 0
z = 0
β = 0.166667

k1 = (d[η] v[η] η (1 - z d[η]) (v[η] - η) - 2 p[η] η (1 - z d[η]) - ϕ[η]^2 d[η] 
     (1 - z d[η]) - Γ p[η] v[η])/((Γ p[η] - (v[η] - η)^2 d[η] (1 - z d[η])) η)

k2 = (d[η] (1 - z d[η]) (v[η] d[η] (v[η] - 2 η) (v[η] - η) + 2 p[η] η + ϕ[η]^2 
     d[η]))/((Γ p[η] - (v[η] - η)^2 d[η] (1 - z d[η])) (v[η] - η) η)

k3 = (p[η] d[η] (2 η (v[η] - η)^2 (1 - z d[η]) + Γ v[η] (v[η] - 2 η) (v[η] - η) +
     ϕ[η]^2 Γ))/((Γ p[η] - (v[η] - η)^2 d[η] (1 - z d[η])) (v[η] - η) η)

k4 = -((ϕ[η] (v[η] + η))/(η (v[η] - η)))

k5 = -(w[η]/(η (v[η] - η)))

sol = NDSolve[{v'[η] == k1, d'[η] == k2, p'[η] == k3, ϕ'[η] == k4, w'[η] == k5, 
   v[1] == (1 - β), d[1] == 1/β, p [1] == (1 - β), ϕ[1] == 0.01, w[1] == 0.02}, 
  {v, d, p, ϕ, w}, {η, 0, 1}, MaxSteps -> 30000]

Please guide me how can I solve the above problem with 4th-order Runge-Kutta method, thanks.

code for RK4 method are given in Solving a system of ODEs with the Runge-Kutta method

but how can I apply those codes to my problem...please guide me...

  • $\begingroup$ Can you show the equations themselves you are trying to solve with the initial/boundary conditions? If you do not know latex, you can scan them from the book and paste the image. $\endgroup$
    – Nasser
    Commented Mar 28, 2015 at 5:53
  • $\begingroup$ @ Nasser... my problem is $${v'[\eta] = k1, d'[\eta] =k2, p'[\eta] = k3, \phi'[\eta] = k4, w'[\eta] = k5, v[1] = (1 - \beta), d[1] = 1/\beta, p [1] = (1 - \beta), \phi[1] = 0.01, w[1] = 0.02}$$ $\endgroup$
    – user45799
    Commented Mar 28, 2015 at 5:55
  • 2
    $\begingroup$ This post has code for Runge-Kutta solving-a-system-of-odes-with-the-runge-kutta-method $\endgroup$
    – Nasser
    Commented Mar 28, 2015 at 6:23
  • 1
    $\begingroup$ After pasting your code here, select them and press Ctrl+K, or add four whitespaces before each line of code. Things like \[Eta] isn't unreadable for Mathematica user, when pasting back to the notebook, it'll get back to the normal greek letters, and if you want to make your code look better in this site (of course it's encouraged), use the link I gave above to convert those "unreadable" symbol to greek letters! $\endgroup$
    – xzczd
    Commented Mar 28, 2015 at 8:35
  • 1
    $\begingroup$ If RunnyKine doesn't appear in the comment under this question, then the "@" won't work. And why are you still adding whitespace between "@" and the name? $\endgroup$
    – xzczd
    Commented Mar 28, 2015 at 9:59

1 Answer 1


According to your statement, I think what you need is just 4th-order Runge-Kutta method, and a completely self-made implementation of 4th-order Runge-Kutta method isn't necessary, then the answer from J.M. has shown you the optimal direction:

(* Unchanged part omitted. *)

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

sol = NDSolve[{v'[η] == k1, d'[η] == k2, p'[η] == k3, ϕ'[η] == k4, w'[η] == k5, 
   v[1] == (1 - β), d[1] == 1/β, p[1] == (1 - β), ϕ[1] == 0.01, w[1] == 0.02}, 
   {v, d, p, ϕ, w}, {η, 0, 1}, 
  Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
    "Coefficients" -> ClassicalRungeKuttaCoefficients}, StartingStepSize -> 1/10000]

However, what I really want to point out is, despite the above code seems to solve your ODE set up to η = 0.0001, I'm afraid it's not reliable at all:

 {{nl, nr}} = (v /. sol)[[1]]["Domain"];
 Plot[{v@n, d@n, p@n, ϕ@n, w@n} /. sol // Evaluate, {n, nl, nr}]

enter image description here

NDSolve by default setting doesn't manage to solve this set of equation, too. It stopped at about η = 0.9576. (I'm not sure what do you mean by saying you have solved it by NDSolve.) I'm not surprised though, your ODEs are non-linear. As for how to solve the ODEs, it's another question. I vote to close this question as a duplicate.

  • $\begingroup$ @xzczd...if answer is not reliable ...please tell me how should I solve my problem....I want to solve this problem numerically for my project... $\endgroup$
    – user45799
    Commented Mar 28, 2015 at 13:20
  • $\begingroup$ @xzczd...when I try to solve the above problem by 'Ndsolve'...is there any command to find which method mathematica has used to solve it $\endgroup$
    – user45799
    Commented Mar 28, 2015 at 13:21
  • 1
    $\begingroup$ @Prasanta Since I'm the author of this post, you don't need to add "@xzczd" to remind me. In my view the most vital part of your problem is, you haven't even find a way to solve your nonlinear ODE. It's easy to create nonlinear ODE, but solving it can be really hard. I suggest you to first make sure if your equations are correct (incorrect translation for the real problem isn't rare!), and consider carefully if you really need your equations to be so complicated. Simplify your model if possible. Diving into the hell of looking for proper combination of options of NDSolve is the last choice. $\endgroup$
    – xzczd
    Commented Mar 29, 2015 at 9:10
  • $\begingroup$ range of η can be between o to 1....as NDSolve stopped at about η = 0.9576. while solving the system.....is there any command to find when we solve the above system by RK4 where it stopped solving .i.e. about η = ? $\endgroup$
    – user45799
    Commented Mar 29, 2015 at 9:39
  • 1
    $\begingroup$ @Prasanta If I correctly understood those material I found, RK4 i.e. classical RK doesn't have an error estimator. Notice the document also gives a example for a fourth order RK called Fehlberg method, which owns an error estimator and can find the singularity at 0.9576. $\endgroup$
    – xzczd
    Commented Apr 1, 2015 at 5:43

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