The vast collection of ODE solvers (e.g., in Odespy) cannot be applied to $$ u^{\prime\prime} + \omega^2 u = 0$$ unless we write this higher-order ODE as a system of 1st-order ODEs.
Introduce an auxiliary variable \( v=u^{\prime} \): $$ \begin{align} u^{\prime} &= v, \tag{1}\\ v^{\prime} &= -\omega^2 u \tag{2} \tp \end{align} $$
Initial conditions: \( u(0)=I \) and \( v(0)=0 \).