function dx=torusDE(t,x)
%Differential equations for the geodesic on a torus
dx=zeros(4,1);

dx(1)=x(2);
dx(2)=-(2+cos(x(1)))*sin(x(1))*x(4)^2;
dx(3)=x(4);
dx(4)=2*((sin(x(1))/(2+cos(x(1)))))*x(2)*x(4);