function dY = motion(theta,Y)

global theta_low theta_high gamma

mu = Y(1);
v  = Y(2);
f  = density(theta); 

dmu = f*(gamma-1/v);
dv  = ((gamma*f)/(gamma*f*theta-2*mu))^2*theta^3;

dY = [dmu ; dv];
