Vectorized Analytic Two Body Propagator (Kepler Universal Variables)

The Keplerian universal variable two body propagation method, using the Newton-Raphson iteration technique, provides a highly efficient propagation routine which doesn't have to transform coordinate systems from cartesian to Keplerian and back again. This method will work on all orbit types (circular, eccentric, parabolic, and hyperbolic) while providing a positional error which is less than 0.5 km (this is obtained from propagation comparisons with the STK implementation of the analytic two body propagation routine).

Example function call is shown below:
>> [r,v] = keplerUniversal(r0,v0,t,mu);

Where:

Inputs:
r0 is the initial position vectors in the ECI coordinate frame of reference (km) [3xN]
v0 is the initial velocity vectors in ECI coordinate frame of reference (km/s) [3xN]
t is the time vector in seconds [1xN]
mu is the Gravitational Constant of the central body which defaults to Earth if not specified (double)

Outputs:
r is the propagated position vector in the ECI coordinate frame of reference (km) [3xN]
v is the propagated velocity vector in ECI coordinate frame of reference (km/s) [3xN]