Solutions to desolve_odeint seem large

edit

I'm trying to figure out what, if anything, is going wrong with this code here. It says in the SageMath documentation that for desolve_odeint, the order of dvars must be the same as that of des. In my case, that means dp/dt=dp, dth/dt=dth, d(dp)/t=g, and d(dth)/t=h. But it's spitting out values that are much larger than what I would expect for t-values between 0 and 10.

t,p,th,dp,dth=var('t,p,th,dp,dth')
g=5*p*dth^2+4*t*dth^2*dp
h=-4*t/(5*p^2+4*t^2)*dth-5*p/(5*p^2+4*t^2)*dp*dth+4*t*dth^3
f=[dp,dth,g,h]
sol=desolve_odeint(f,[5,0,1,1],srange(0,10),[p,th,dp,dth],ivar=t)
sol[2]