//A rossler iterator
#include <stdlib.h>
#include <math.h>
#include <fstream.h>
#include <iostream.h>
#include "strstream.h"
#include "nr.h"
#include "nrutil.h"
#include "nrutil.c"
//#include "rk4.c"

void derivsy(double, double[], double[]);

int main(int argc, char *argv[]){
  int i1,i2,i3,j,n1,n,N;  double x,a,b,r,fixedx,fixedy,fixedz; double h=0.001; double *y; 
 double *dydx; double *yout;
 n=3;
 //a=0.2;
 a=10; b=((1.0)*8)/3; 
 char *str = (char*)malloc(40);
 

 for(r=28; r<60; r+=1){

 sprintf(str,"%f",r);



 y = dvector(1,n);  dydx = dvector(1,n); yout = dvector(1,n);


 //Initial values
 y[1]=0.25; y[2]=0.0; y[3]=0.1;


//this first main loop is here to remove original transients *
 for(i2=0; i2 <=10000; i2++){

 dydx[1] = a*(y[2]-y[1]); //-y[2]-y[3];
 dydx[2] = r*y[1] - y[2] -y[1]*y[3]; //y[1] + a*y[2];
 dydx[3] = y[1]*y[2] - b * y[3];

   y[1]+=h*dydx[1];
   y[2]+=h*dydx[2];
   y[3]+=h*dydx[3];}

 ofstream outf(str);

 for(i2=1; i2<=200000; i2++){ 

  dydx[1] = a*(y[2]-y[1]); //-y[2]-y[3];
 dydx[2] = r*y[1] - y[2] -y[1]*y[3]; //y[1] + a*y[2];
 dydx[3] = y[1]*y[2] - b * y[3];


    y[1]+=h*dydx[1];
    y[2]+=h*dydx[2];
    y[3]+=h*dydx[3];

    outf<<" "<<y[1]<<" "<<y[2]<<" "<<y[3]<<"\n";
 }  //ends i2 loop

 fixedx=sqrt(b*(r-1));
 fixedy=sqrt(b*(r-1));
 fixedz=r-1;

 //Print out fixed points;
 outf<<" "<<fixedx<<" "<<fixedy<<" "<<fixedz<<"\n";
 outf<<" "<<-fixedx<<" "<<-fixedy<<" "<<fixedz<<"\n";
 //outf<<" "<<b<<" "<<(r-1)<<" "<<sqrt(b*(r-1))<<" "<<"\n";
 outf.close();

 } //ends a;

}  //ends main