/* file model.c  4.7.1995 

function rhs 
evaluates the right-hand-side functions 
    f(x) : Rdim -> Rdim
of the equation
    x' = f(x)

this demo is written to compute the Lyapunov exponents of the Lorenz system

*/

rhs (time, vec, pars, funs, dummy1, dummy2, dim)
 real time;
 /* independent variable (used for non-autonomous systems only), input */
 vector
  vec,      /* vector of state variables, input  */
  pars,     /* vector of parameters,      input  */
  funs,     /* vector of rhs functions,   output */
  dummy1 ;  /* not used here, used for stiff systems */
 matrix dummy2;  /*  dtto */

 int dim;   /* dimension of the full system 
                e.g. for integrating 3 Lorenz equations ... dim = 3
                         3D Lorenz with 1 L.E.                    6
                                        2 L.E.                    9
                                        3 L.E.                   12
*/

{
 real 
  x = vec [0],
  y = vec [1],
  z = vec [2],

  s = pars[0],
  r = pars[1],
  b = pars[2];

 matrix J; /* matrix of partial derivatives of funs  d funs  / d vec */

 funs [0] = s*(y-x);
 funs [1] = r*x-y-x*z;
 funs [2] = x*y-b*z;

 if (dim > 3) {
  J[0][0] = -s;
  J[0][1] =  s;
  J[0][2] =  0;
 
  J[1][0] =  r - z;
  J[1][1] = -1;
  J[1][2] = -x;
 
  J[2][0] =  y;
  J[2][1] =  x;
  J[2][2] = -b;
 };

if (dim > 3) timesmv (funs+3,J,vec+3,3,3);   
if (dim > 6) timesmv (funs+6,J,vec+6,3,3);   
if (dim > 9) timesmv (funs+9,J,vec+9,3,3);   

} /* rhs_Lor */
/*********************************************/