# RK4 of second order ODEs system

I'm new here and I have to program explicitly a Runge-Kutta 4 Method for a second order ODEs $$10 \times 10$$ system: $$X''(t)=AX(t)+F(t)$$, and the instruction was to make a change of variable $$y(t)=x'(t)$$, expanding it into a $$20 \times 20$$ first order system: $$Y'(t)=BY(t)+G(t)$$. The $$10 \times 10$$ system works good when I solved it using NDSolve, and also I have already programed RK4 for a particular case of a first order ODE, in this code $$x'(t)=-\frac{1}{2}x(t)-x^2(t),\;x(t=1)=0.1$$:

t0 = 1;
x0 = 0.1;
tf = 30;
n = 100;
h = (tf - t0)/n;
x[0] = x0;
t[n] = tf;
f[t_, x_] = -1/2 x/t - x^2;

For[i = 0, i < n, i++,
{
t[i] = t0 + i*h;
k1 = f[t[i], x[i]];
k2 = f[t[i] + (1/2)*h, x[i] + (1/2)*h*k1];
k3 = f[t[i] + (1/2)*h, x[i] + (1/2)*h*k2];
k4 = f[t[i] + h, x[i] + h*k3];
x[i + 1] = x[i] + (1/6)*h*(k1 + 2*k2 + 2*k3 + k4);
}
];
A = ListPlot[{Table[{t[i], x[i]}, {i, 1, n} ] }, PlotLegends -> {RK4}]


And I don't know how to use it with my code for the 20x20 system:

omega = 3;
B = {{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {-4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {2, 0, -4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 2, 0, -4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 2, 0, -4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 2, 0, 4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 4, 0, 2, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 4, 0, 2, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 4, 0, 2, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 4, 0, 2,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, -2, 0}};
Y[t_] = {x1[t], y1[t], x2[t], y2[t], x3[t], y3[t], x4[t], y4[t],
x5[t], y5[t], x6[t], y6[t], x7[t], y7[t], x8[t], y8[t],
x9[t], y9[t], x10[t], y10[t]};
G[t_] = {0, Cos[omega*t], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
Y[0]=={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}
Sys = B.Y[t] + G[t]


I'd appreciate if someone helps me, I've tried in many ways and I couldn't do it.

P.S.: I'm Chilean and my English is not very good, so sorry if something is misspelled.

• 1. It's not a $20\times 20$ system, $B$ is indeed a $20\times 20$ matrix though. 2. Why is there a x/t in the definition of f? 3. "I've tried in many ways and I couldn't do it." Please show us how you tried so we can give specific suggestion. – xzczd Oct 5 at 8:50
• If we write the first-order system in the form $Y' = f(t, Y)$, where $Y$ is a 20-vector and $f$ returns a 20-vector, then you should be able to use your For loop. Then $f$ should be f[t_, y_] := B.y + G[t]. – Michael E2 Oct 5 at 15:02