The 2nd ODE that I want to solve numerically is:

xddot = ( 0.5 * rho * A * Cl * ( xdot^2 + vy^2) ) / m

ODE

The vy term is a sinusoidal function and does not need to be solved numerically.

I have transfromed it into a coupled 1st ODE by:

xdot = vx

vxdot = ( 0.5 * rho * A * Cl * (vx^2 + vy^2) ) / m

The fact that the xdot term is squared in the ODE makes the solution very unstable. What are possible solutions to mitigate this?

Attached is the output of the ODE that I am getting.

Output