function dYdt = kepler_ode(t,Y)
% assume 'Y' is a vector of form y = [ x ; y ; vx ; vy ]
r = [Y(1);Y(2)];
v = [Y(3);Y(4)];
% use the following three variables from the main worksheet
% must declare them 'global' here and there
global mStar mPlanet G
% assumes the sun is located at r=(0 0)
F = - mStar * mPlanet * G / norm(r)^3 .* r;
a = F ./ mPlanet;
dYdt = [ v(1) ; v(2) ; a(1) ; a(2)]; % must column vector