JacobianODE:

--------------------------------------------------------------------------
   This function computes the Jacobian matrix for a right hand-side.
   The central difference approximation good to O(h^2) is used.
   For any function f(x,t) it will compute the f0 and a matrices where
   dx/dt = ax + f0.

   f(x,t) ~= f(xOP,tOP) + a(xOP,tOP) x + ...

   funFcn is input as a character string 'xxxxx' and must be of the form
   xdot = xxxxx ( t, x, flag, d )
   which is the same form used in ode113. 

--------------------------------------------------------------------------
   Form:
   [a, fOP] = JacobianODE( funFcn, tOP, xOP, d )
--------------------------------------------------------------------------

   ------
   Inputs
   ------
   funFcn         (1,:)  T'function name'
   tOP            (1,1)  Time
   xOP            (n,1)  State at the operating point
   d                     Optional data structure

   -------
   Outputs
   -------
   a              (n,n)  Jacobian matrix of first partials       
   fOP            (n,1)  f(xOP,tOP) at the operating point

--------------------------------------------------------------------------