function f = satdyn( t,state )

xp = state(1);  % x position body 
yp = state(2);  % y position body 
xv = state(3);  % x velocity body 
yv = state(4);  % y velocity body 
mu = 4e14;
r3 = (xp.^2 + yp.^2) * sqrt(xp.^2 + yp.^2);
xa = -mu/r3*xp;
ya = -mu/r3*yp;
f(1) = state(3);
f(2) = state(4);
f(3) = xa;
f(4) = ya;
f=f';