How to make Runge-Kutta for solving nonlinear ODE system in Mathematica

Question

How to make Runge Kutta method for system of non linear equations in this Matrix form. Matrices A3, B3 and B4 are functions of matrix X3. Initial conditions are X1=X2=X3=X4=0 and system is here

$$ \frac{d}{dt}\left( \begin{array}{c} \text{X1} \\ \text{X2} \\ \text{X3} \\ \text{X4} \end{array} \right)=\left( \begin{array}{cccc} -\text{A2} & -\text{A1} & -\text{A3} & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ -\text{B3} & 0 & -\text{B2}-\text{B4} & \text{B1} \end{array} \right).\left( \begin{array}{c} \text{X1} \\ \text{X2} \\ \text{X3} \\ \text{X4} \end{array} \right)+\left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \text{B5} \end{array} \right).\left( \begin{array}{cccc} 0 & 0 & 0 & 1 \end{array} \right) $$

$$ \text{X1}=\left( \begin{array}{c} \text{x11} \\ \text{x12} \\ \text{x13} \end{array} \right); $$

$$ \text{X2}=\left( \begin{array}{c} \text{x21} \\ \text{x22} \\ \text{x23} \end{array} \right); $$

$$ \text{X3}=\left( \begin{array}{c} \text{x31} \\ \text{x32} \\ \text{x33} \end{array} \right); $$

$$ \text{X4}=\left( \begin{array}{c} \text{x41} \\ \text{x42} \\ \text{x43} \end{array} \right); $$

Where we have matrices

$$\text{A1}=\left( \begin{array}{ccc} \frac{640576501799}{76356} & 0 & \frac{311934265235}{57839} \\ 0 & \frac{1285802795705}{40871} & 0 \\ \frac{388065734765}{30838} & 0 & \frac{523620702496}{6935} \end{array} \right);$$

$$ \text{A2}=\left( \begin{array}{ccc} \frac{980000000000000}{1168149} & 0 & \frac{210000000000000}{389383} \\ 0 & \frac{1225000000000000}{389383} & 0 \\ \frac{490000000000000}{389383} & 0 & \frac{2940000000000000}{389383} \end{array} \right);$$

$$ \text{A3}=\left( \begin{array}{ccc} 0 & \frac{152600000000000000 \text{x31}}{14532941709}+\frac{76300000000000000 \text{x33}}{4844313903} & 0 \\ 0 & \frac{953750000000000000 \text{x32}}{14532941709} & 0 \\ 0 & \frac{76300000000000000 \text{x31}}{4844313903}+\frac{1068200000000000000 \text{x33}}{4844313903} & 0 \end{array} \right);$$

$$ \text{B1}=\left( \begin{array}{ccc} \frac{5874951206}{12955399} & 0 & \frac{3212738087}{4414317} \\ 0 & \frac{8215070163}{2795011} & 0 \\ \frac{8400845503}{6412667} & 0 & \frac{6246501886}{520165} \end{array} \right); $$

$$ \text{B2}=\left( \begin{array}{ccc} \frac{14850000000000}{327471103} & 0 & \frac{71500000000000}{982413309} \\ 0 & \frac{96250000000000}{327471103} & 0 \\ \frac{42900000000000}{327471103} & 0 & \frac{393250000000000}{327471103} \end{array} \right); $$

$$ \text{B3}=\left( \begin{array}{ccc} 0 & -\frac{180000000000000000 \text{x31}}{1027581737}-\frac{40000000000000000 \text{x33}}{79044749} & 0 \\ 0 & -\frac{700000000000000000 \text{x32}}{440392173} & 0 \\ 0 & -\frac{40000000000000000 \text{x31}}{79044749}-\frac{660000000000000000 \text{x33}}{79044749} & 0 \end{array} \right); $$

$$ \text{B4}=\left( \begin{array}{ccc} \frac{55687500000000000000 \left(\frac{16 \text{x31}^2}{11781}-\frac{640 \text{x31} \text{x33}}{2909907}\right)}{327471103}+\frac{160875000000000000000 \left(\frac{16 \text{x31} \text{x33}}{11781}-\frac{640 \text{x33}^2}{2909907}\right)}{327471103} & \frac{540900000000000000000 \text{x31} \text{x32}}{9625358130479}+\frac{120200000000000000000 \text{x32} \text{x33}}{740412163883} & \frac{55687500000000000000 \left(-\frac{640 \text{x31}^2}{2909907}+\frac{2456 \text{x31} \text{x33}}{14549535}\right)}{327471103}+\frac{160875000000000000000 \left(-\frac{640 \text{x31} \text{x33}}{2909907}+\frac{2456 \text{x33}^2}{14549535}\right)}{327471103} \\ \frac{505312500000000000000 \left(\frac{16 \text{x31} \text{x32}}{11781}-\frac{640 \text{x32} \text{x33}}{2909907}\right)}{327471103} & \frac{2103500000000000000000 \text{x32}^2}{4125153484491} & \frac{505312500000000000000 \left(-\frac{640 \text{x31} \text{x32}}{2909907}+\frac{2456 \text{x32} \text{x33}}{14549535}\right)}{327471103} \\ \frac{160875000000000000000 \left(\frac{16 \text{x31}^2}{11781}-\frac{640 \text{x31} \text{x33}}{2909907}\right)}{327471103}+\frac{2654437500000000000000 \left(\frac{16 \text{x31} \text{x33}}{11781}-\frac{640 \text{x33}^2}{2909907}\right)}{327471103} & \frac{120200000000000000000 \text{x31} \text{x32}}{740412163883}+\frac{1983300000000000000000 \text{x32} \text{x33}}{740412163883} & \frac{160875000000000000000 \left(-\frac{640 \text{x31}^2}{2909907}+\frac{2456 \text{x31} \text{x33}}{14549535}\right)}{327471103}+\frac{2654437500000000000000 \left(-\frac{640 \text{x31} \text{x33}}{2909907}+\frac{2456 \text{x33}^2}{14549535}\right)}{327471103} \end{array} \right); $$

$$ \text{B5}=\left( \begin{array}{c} \frac{37125000}{3241} \\ 0 \\ \frac{107250000}{3241} \end{array} \right); $$

Answer 1

Function NDSolve in Mathematica has the Runge-Kutta methods built-in:

NDSolve[{x1'[t] == -x1[t] x2[t]^2, x2'[t] == -x1[t] (1 + x2[t]^2), 
      x1[0] == 1, x2[0] == 1}, {x1, x2}, {t, 0, 2}, 
      Method -> "ExplicitRungeKutta"]

According to NDSolve's online reference page, method "ExplicitRungeKutta" gives explicit Runge-Kutta methods with adaptive embedded pairs of 2(1) through 9(8), method "ImplicitRungeKutta" gives families of arbitrary-order implicit Runge-Kutta methods.

Order of the Runge-Kutta method can be set as follows:

NDSolve[{x1'[t] == -x1[t] x2[t]^2, x2'[t] == -x1[t] (1 + x2[t]^2), 
  x1[0] == 1, x2[0] == 1}, {x1, x2}, {t, 0, 2}, 
 Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 7}]

NDSolve permits writing your own method, and an example of implementing classic Runge-Kutta method of order 4 is given in online tutorial "NDSolve Method Plugin Framework".

---

Added: This example solve a vector system in NDSolve using Runge-Kutta method:

xvec = {x1[t], x2[t], x3[t]};
x0vec = {1, 2, 3};
amat = {{0, -x3[t], x2[t]}, {-x3[t], 0, x1[t]}, {1, 1, 1}};
NDSolve @@ {Flatten@{Thread[D[xvec, t] == amat.xvec], 
    Thread[x0vec == (xvec /. t -> 0)]}, {x1, x2, x3}, {t, 0, 2}, 
  Method -> "ExplicitRungeKutta"}