% This code is copied from p. 201 of "Inverse Problems" by C.W. Groetsch (MAA, 1999).
%
% The function defined below sets up the equations of the motion in the inverse-square field:
%
%   d^2 r / dt^2 = - r /|r|^3,   where r=<x,y>.
%
%  This equation is rewritten as a system of first-order differential equations
%  for the components of the  "vector"  z=[x,y,xdot,ydot].
%
% This function is invoked by  hw2_p7.m (Spring 2013). 

function  zdot=hw2_p7_funorbit(t,z)

zdot(1) = z(3);           % dot(x) = xdot
zdot(2) = z(4);           % dot(y) = ydot
zdot(3) = -z(1)/sqrt(z(1)^2+z(2)^2)^3;   % dot(dot(x)) = -x/(x^2+y^2)^(3/2)
zdot(4) = -z(2)/sqrt(z(1)^2+z(2)^2)^3;   % dot(dot(y)) = -y/(x^2+y^2)^(3/2)
zdot=zdot';       % this makes the row-vector zdot defined above into a column-vector
