Hello, @F. Semih Dündar! I can see a couple of problems with your code. First of all, there is no need to define $q$ and $qdot$ as formal functions, since that is implicitly assumed. Second, you only need to pass the RHS (right hand side) of the differential equations. Finally, the `ics`

parameter must be of the form `[t0, q0, qdot0]`

.

With those considerations, your code should look like this:

```
var('t, qdot, q') # independent variable and functions
# Equations of motion
eom1 = qdot
eom2 = -(qdot(t))^2 - q(t)
# Initial conditions
q0 = 0
qdot0 = 1
sol = desolve_system_rk4([eom1, eom2], ivar=t, vars=[q, qdot], ics = [0, q0, qdot0], end_points=1.8, step=0.01)
S = [(t, q) for t, q, qdot in sol]
list_plot(S, plotjoined=True)
```

I think a few additional clarifications are in order. As mentioned, you only need the RHS of the DEs; that is why I removed part of $oem1$ and $oem2$ (specifically, the LHS and the `==`

). The `vars`

parameter implicitly indicate that your system is of the form

$$q' = oem1$$
$$qdot' = oem2$$

so this code defines the same system you intended.

I have used the `end_points=1.8`

parameter in order to solve in the interval $[0, 1.8]$---after $t\approx1.8$, the value of $q$ is negative, which I assumed not to be desirable, although I don't know the context in which you are trying to solve this (so maybe you would like to modify that parameter if needed).

On the other hand, when I plotted the solution, due to the small domain, I obtained a pointy curve, so I used the `step=0.01`

parameter, in order to decrease the step length, and obtain a more refined solution.

The last two lines are intended to obtain a plot of $q$ alone, without including $qdot$. (I assume they are clear to you, but feel free to ask otherwise.)

Here is the resulting image:

I hope this helps!