

function F=galactic_ODEs(t,Y,flag)

global k m L0 L g

x=Y(1);
y=Y(2);
xp=Y(3);
yp=Y(4);

sq = sqrt((L-y).^2+x.^2);
L0msqosq = (L0-sq)./sq;

xpp = (k/m)*x.*L0msqosq;
ypp = -g+(k/m)*(y-L).*L0msqosq;

F=[xp;yp;xpp;ypp];




