EOM1JAC.m

function dA_dt = EOM1JAC(t, A, varargin)
  global indexes

  r  = A(indexes.r);
  v  = A(indexes.v);
  J  = reshape(A(7:42,end), 6, 6);

  r2 = dot(r,r);
  rm = sqrt(r2);
  r3 = r2 * rm;
  r5 = r2 * r3;

  % propagate the states

  rdot  = v;
  vdot  = - r / r3;

  dX_dt = [ rdot' vdot' ]';

  % propagate the jacobian

  dR_dR = zeros(3,3);
  dR_dV = eye(3,3);
  dV_dR = 3/r5*r*r' - eye(3,3)/r3;
  dV_dV = zeros(3,3);

  dJ_dt = [ dR_dR dR_dV; dV_dR dV_dV ] * J;

  % assemble the vector of differentials

  dA_dt = [
    dX_dt
    dJ_dt(:) % ordered by columns
    ];
end