constant_control_action.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# constant control action:
u = 0.5


# define the system dynamics
def system_dynamics(_t, z):
    return [z[1],
            -z[0]*z[1] + u]


t_final = 4
z_init = [1, 1]
solution = solve_ivp(system_dynamics, [0, t_final], z_init,
                     t_eval=np.linspace(0, t_final, 100))

plt.plot(solution.t, solution.y.T)
plt.show()