% Equations of motion in Earth's 
% gravity well (Earth at origin).
% q is a column vector, [x; y; vx; vy], 
% in mks units.
function dqdt = earth(q, t)

g = 9.81;   % Earth surface gravity (m/s)
R = 6371e3; % Earth radius (m)

r = q(1:2); % position vector
v = q(3:4); % velocity vector

% Equations of motion
drdt = v;
dvdt = - g * r * (R^2) / norm(r)^3;

dqdt = [drdt; dvdt];