#include <math.h>

#include "chapter8.cpp"

typedef double runtype;
ofstream outfile1("euler.txt");
ofstream outfile2("taylor.txt");
ofstream outfile3("rk4.txt");
ofstream outfile4("rk45.txt");
ofstream outfile5("rk45ad.txt");

inline runtype f1(runtype t, runtype x){ return(1 + x*x + t*t*t); }
inline runtype f2(runtype t, runtype x) { return 2 + (x - t - 1)*(x - t - 1); }
inline runtype f3(runtype t, runtype x) { return 2.0 + (x - t - 1.0)*(x - t - 1.0); }
inline int sign(runtype x) { return (x >= 0.0) ? 1 : -1; }
inline runtype max(runtype x, runtype y) { return (x >= y) ? x : y; }
inline runtype f4(runtype t, runtype x) { return 3 + 5*sin(t) + 0.2*x; }


void main()
{
  int n = 100;
  runtype a = 1.0, b = 2.0, x = -4.0;
  int k;

  euler(a, b, x, n);

  a = 1.0, b = 2.0, x = -4.0;
  taylor(a, b, x, n);

  n = 72;
  a = 1.0, b = 1.5625, x = 2.0;
 
  runtype h, t;
  h= (b-a)/n;
  t = a;
  rk4(outfile3, f2, t, x, h, n);

  runtype epsilon;

  outfile4 << "k=  0:x(" << setw(10) << t << ") = " << setw(10) << x << endl;
  for (k = 1; k <= n;  k++)
  {
    rk45(f3, t, x, h, epsilon);
	outfile4 << "k =" << setw(3) << k << ": x(" << setw(10) << t <<") =" << setw(10) << x 
		<<"   e =" << setw(10) << epsilon<< endl;
  }

  const int itmax = 1000;
  const runtype epsmin = 0.00000001, 
              epsmax = 0.00001, 
              hmin = 0.000001, 
              hmax = 1.0;
  runtype tb; 

  t = 0.0;
  x = 0.0;
  h = 0.01;
  tb = 10.0;
  rk45ad (outfile5, f4, t, x, h, tb, itmax, epsmin, epsmax, hmin, hmax);
}