function [Xdot] = twobodyEOM(t,X)
% =========================================================================

% This is your gravitational parameter of your center body 
% mu = GM of Earth 398,600 km^3/s^2
mu = 398600;

% This is your states: x, y, z, xdot, ydot, zdot of your s/c
x = X(1);  
y = X(2);  
z = X(3);
xdot = X(4);
ydot = X(5);
zdot = X(6);

% This is your radius 
r = sqrt( x^2 + y^2 + z^2 );  
% An easier and less computationaly intensive way is to write r = norm(X(1:3))

% Here are your equations of motion
accelx = - mu * x / r^3;
accely = - mu * y / r^3;
accelz = - mu * z / r^3;

% Here are your outputs of all the velocity and acceleration
Xdot    = zeros(6,1);
Xdot(1) = xdot;
Xdot(2) = ydot;
Xdot(3) = zdot;
Xdot(4) = accelx;
Xdot(5) = accely;
Xdot(6) = accelz;

% =========================================================================