// LorenzXZ.c: Lorenz simulator to be run under ch.

#include <stdio.h>
#include <math.h>
#include <numeric.h>

#define NVAR 3
#define POINTS 50000

void lorenz(double t, double y[], double dydt[]);

int main()
  {
  double t0=0, tf=100, y0[NVAR]={1,0,0};
  double t[POINTS], y[NVAR][POINTS];
  int i, points;
  /* return -1: divergent return >1: ok, num of points */
  points = odesolve(t, y, lorenz, t0, tf, y0);
  if(points > 0)
	plotxy(y[0], y[2], "Lorenz system", "x", "z");
  else
	printf("odesolve() failed\n");
  }

void lorenz(double t, double y[], double dydt[])
  {
  double sigma=10, b=8/3, r=24.5;
  dydt[0] = sigma*(y[1]-y[0]);
  dydt[1] = r*y[0] - y[1] -y[0]*y[2];
  dydt[2] = y[0]*y[1] - b*y[2];
  }