### Solving 2nd order ode with desolve_system_rk4

Dear All,

I have the following nonlinear ode to solve numerically:

$$\ddot q(t) = -\dot q(t)^2 - q(t)$$.

Since desolve_system_rk4 cannot solve 2nd order ode I make this ode a system of two 1st order odes in the following way:

$$q'(t) = qdot(t)$$
$$qdot'(t) = -qdot(t)^2 - q(t)$$

However desolve_system_rk4 cannot still solve this 2 odes numerically. It raises an error. Any help is appreciated.

Here is the code I wrote:

```
t = var('t')
q = function('q')(t)
qdot = function('qdot')(t)
# Equations of motion
eom1 = q.diff(t) == qdot(t)
eom2 = qdot.diff(t) == -(qdot(t))^2 - q(t)
# Initial conditions
q0 = 0
qdot0 = 1
sol = desolve_system_rk4([eom1, eom2], [q, qdot], ivar = t , ics = [q0, qdot0])
```