I'm trying to find lypunov exponent for this  system of ODEs. I know I need to take the Jacobian but not sure if it's possible the way it's currently defined. If anyone could provide insight it would be much appreciated.

Eqns:= diff(omega(t),t)=-(G*MSat*beta^(2)*(xH(t)*sin(theta(t))-yH(t)* cos(theta(t)))*(xH(t)*cos(theta(t))+yH(t)*sin(theta(t))))/((xH(t)^(2)+yH(t)^(2))^(2.5)),diff(theta(t),t)=omega(t), diff(xH(t),t)=vxH(t),diff(vxH(t),t)=-(G*M*xH(t))/((xH(t)^(2)+yH(t)^(2))^(1.5)),diff(yH(t),t)=vyH(t),diff(vyH(t),t)=-(G*M*yH(t))/((xH(t)^(2)+yH(t)^(2))^(1.5)): ;

ICs := omega(0) = omega0, theta(0) = theta0, xH(0) = a*(1+e), yH(0) = 0, vxH(0) = 0, vyH(0) = sqrt(G*M*(1-e)/(a*(1+e)));

I want the exponent for omega, I  procedure that takes some initial conditions, changing just w0, and computes the long term value of omega. This plots a sort of bifurcation diagram. I'd like an estimate of the exponent to compare what I see. 

thanks for the help

ft