function yp = f(t,y)

% y= x1 x2 x3 x4 l1 l2 l3 l4 

global mass momin 

momineff=momin+mass*y(2)^2;

lam1=y(5);
lam2=y(6);
lam3=y(7);
lam4=y(8);

u1=lam3/momineff;
u2=lam4/mass;

yp(1)=y(3);
yp(2)=y(4);
yp(3)=(u1-2*mass*y(2)*y(3)*y(4))/momineff;
yp(4)=u2/mass+y(2)*y(3)^2;

yp(5)=0;
yp(6)=2*mass/momineff*(y(3)*y(4)+y(2)*yp(3))*lam3-y(3)^2*lam4;
yp(7)=-lam1+2*mass/momineff*y(2)*y(4)*lam3-2*y(2)*y(3)*lam4;
yp(8)=-lam2+2*mass/momineff*y(2)*y(3)*lam3;

