Systems of second order ODEs can be reduced to systems of first order ODEs.

For example let us start with this system of second order ODEs:

- $y_1'' + y_2 + 1 = 0$
- $y_2'' + y_1 = 0$

Introduce additional functions:

- $y_3 = y_1'$
- $y_4 = y_2'$

We obtain a system of first order ODEs:

- $y_1' = y_3$
- $y_2' = y_4$
- $y_3' = -y_2 - 1$
- $y_4' = -y_1$

Sage numerical solution:

```
sage: T = ode_solver()
sage: f = lambda t, y: [y[2], y[3], -y[1] - 1, -y[0]]
sage: T.function = f
sage: T.ode_solve(y_0=[1, 1, 0, 0], t_span=[0, 20], num_points=1000)
```

Plot $y_1$:

```
sage: f = T.interpolate_solution()
sage: plot(f, 0, 5).show()
```

Plot $y_2$:

```
sage: f = T.interpolate_solution(i=1)
sage: plot(f, 0, 5).show()
```