/* C program that computes the Lyapunov exponents
   for the chaotic chemical system from

J.L.Hudson & O.E.Rossler:
Chaos in Simple Three- and Four-variable Chemical Systems
in Lecture Notes in Biomathematics 55:
Modelling of Patterns in Space and Time
pp.135-145
eds. W.Jager & J.D.Murray
Springer-Verlag Berlin 1984

To run it 
(1) save this code in a file hudson_le.c
(2) save easynum library in a file easynum.c
(3) compile it by
$ cc hudson_le.c -lm
(4) run it by
$ a.out

*/
# include "easynum.c" /* http://staff.vscht.cz/~pokornp/easynum */
/******************************/
void hudson (t,y,p, f,g,J,d)
 real t;
 vector y,p, f,g;
 matrix J;
{
 real 
  k1=0.01, k2=0.02, k3=1.0, k=0.08, k4=0.11, k5=0.0005,

  a=y[0],
  b=y[1],
  c=y[2];

 f[0] = k1-k2*(a-c)-k3*b*a;
 f[1] = k3*a*b-k4*b/(b+k)+k5;
 f[2] = k2*(a-c);

 J[0][0] = -k2 - k3 * b;
 J[0][1] = -k3 * a;
 J[0][2] = k2;
 J[1][0] = k3 * b;
 J[1][1] = k3 * a - k4 * k / ((b+k)*(b+k));
 J[1][2] = 0;
 J[2][0] = k2;
 J[2][1] = 0;
 J[2][2] = -k2;

 if (d>3) timesmv (f+3,J,y+3,3,3);
 if (d>6) timesmv (f+6,J,y+6,3,3);
 if (d>9) timesmv (f+9,J,y+9,3,3);
} /* hudson */
/******************************/
int main ()
{
 real TT=1000, T=100000;
 vector
  x,           /* vector of the state variables */
  le,          /* vector of resulting Lyapunov exponents */
  pars;        /* vector of parameters */
 int
  output_wanted = 0, /* try 1 to see the results of integration  */
  steps = 10,  /* number of integration steps between re-orthonormalizations */
  nle = 3,     /* number of L.E. wanted */
  d = 3;       /* dimension of the original system
                  dimension of the full system will be dim = d * (nle+1)   */
 x[0] = 0.8;
 x[1] = 0.01;
 x[2] = 1.0;

 /* compute the Lyapunov exponents */
 LE (x,le,hudson,pars,T,TT,output_wanted,steps,nle,d);

 /* print them */
 writevector ("LEs are: ",le,nle);
 return (0);
}
/*********************************************/