function yprime = f(tk, yk)
% The system can be written as (rx = x-coordinate of r etc):
% rx' = vx   v = velocity
% ry' = vy
% rz' = vz
% vx' = c sum_j=1^n (rx - px_j) / || r - p_j ||^3
% vy' = c sum_j=1^n (ry - py_j) / || r - p_j ||^3
% vz' = c sum_j=1^n (rz - pz_j) / || r - p_j ||^3
%
% In Matlab vi store rx, ry in the vector r, so rx = r(1) etc.
% Let y = [rx ry rz vx vy vz]' i.e. y = [r; v];

global P c

r = yk(1:3);  % Pick out r and v
v = yk(4:6);

s = [0; 0; 0];
for j = 1:size(P, 2)
  d = r - P(:, j);
  s = s + d / norm(d)^3;
end
yprime = [v; c * s];
