Jacobian:

--------------------------------------------------------------------------
   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.

   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 ( x, t, {p1,p2,p3,p4,p5,p6,p7,p8,p9,p10})
   which is the same form used in RK2 and RK4. tOP is optional. 
   If not needed pass [].

--------------------------------------------------------------------------
   Form:
   [a, fOP] = Jacobian( funFcn, xOP, tOP, varargin )
--------------------------------------------------------------------------

   ------
   Inputs
   ------
   funFcn                'function name' or handle
   xOP            (n,1)  State at the operating point
   tOP                   Time
   varargin              Optional arguments

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

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