Keplerian State Transition Matrix

Given initial positions and velocities, calculate the keplerian state transition matrix for an orbit

Phi = keplerSTM(x0,dt,mu) will return the state transition matrix Phi
for a set of objects described by initial positions and velocities (x0)
with gravitational parameters mu, propogated along Keplerian orbits for
time dt. The positions and velocities at time dt is then given by Phi*x0.

x0 must be a row vector of length 6n where n is the number of planets.
For each planet, x0 must include the 3 position and 3 velocity values
in the order [r1,r2,r3,v1,v2,v3].' These positions and velocities are
measured in a cartesian coordinate system with the central object
(i.e. star) at the origin. For multiple planets, simply stack
several of these vectors on top of each other. mu must contain n
values equal to G(m+ms) where G is the gravitational constant, m is the
mass of the orbiting object, and ms is the mass of the central object.
dt is a scalar value of time to propogate.
NOTE: All units should be consistent. If the positions are in AU and
velocities are in AU/day, then dt must be in days, and mu must be in
AU^3/day^2.

This is calculated using the algorithm described in Shepperd, 1984
which employs Goodyear's universal variables and solves the Kepler
problem using continued fractions.

Example:
%propogate a planet 10 days along its orbit:
x0 = [-0.8684, -0.6637, 0.7207, 0.0039, 0.0056, 0.0120].';
x1 = keplerSTM(x0,10,2.6935e-04)*x0;