/* This program integrates the Lorenz system */ 
# include "easynum.c" /* http://staff.vscht.cz/~pokornp/easynum/ */
/************************************/
void lorenz (time,xyz,pars, xyzdot,g,J,dim)
 real time;
 vector xyz,pars,xyzdot,g;
 matrix J; 
 int dim;
{
 real sigma = 10, b = 8.0/3.0, r = 28, x = xyz[0], y = xyz[1], z = xyz[2]; 

 xyzdot[0] = sigma * (y-x);
 xyzdot[1] = r * x - y - x * z;
 xyzdot[2] = x * y - b * z;
} /* lorenz */
/************************************/
int main ()
{
 real t = 0,           /* time */ 
  Ts=0.01              /* sampling period */;
 vector xyz,pars;      /* state vector and parameters */
 int i,imax=10000,dim=3;
 
 xyz[0] =  6;          /* set initial conditions */
 xyz[1] = 10;
 xyz[2] = 14;

 for (i=0;i<imax;i++) {
  integdt (&t,Ts,xyz,lorenz,pars,dim);  /* integrate lorenz for time Ts */
  printf ("%G %G %G %G\n",t,xyz[0],xyz[1],xyz[2]);
 };
 return (0);
} /* main */
/************************************/