#include "chapter12.cpp"

typedef double runtype;
typedef std::vector<runtype> Array;
typedef std::vector<Array> Matrix;



inline runtype u(runtype x)
{
  return (exp(x) - 3 * sin(x));
}

inline runtype v(runtype x)
{
  return (-1.0);
}

inline runtype w(runtype x)
{
  return (1.0);
}

ofstream outfile1("bvp1.txt");
ofstream outfile2("bvp2.txt");

void main()
{
  const int n = 99;
  const runtype ta = 1.0,tb = 2.0,alpha = 1.09737491,beta = 8.63749661;
  Array a(n+1,0), b(n+1,0), c(n+1,0),d(n+1,0), y(n+1,0);
  int i;
  runtype error, h, t;

// BVP1
  
  h = (tb - ta) / n;
  for (i = 1 ; i < n ; i++)
  {
    t = ta + i * h;
    a[i] = - (1.0 + (h/2) * w(t));
    d[i] = 2.0 + h * h * v(t);
    c[i] = -(1.0 - (h/2) * w(t));
    b[i] = - h * h * u(t);
  }

  b[1] -= a[1] * alpha;
  b[n-1] -= c[n-1] * beta;
  for (i = 1; i < n; i++)
    a[i] = a[i+1];
  tri(n-1, a, d, c, b, y);
  error = exp(ta) - 3.0 * cos(ta) - alpha;
  outfile1 << setw(15) << "ta" << setw(15) << "Solution" << setw(15) << "Error" << endl;
  outfile1 << setw(15) << ta << setw(15) << alpha << setw(15) << error << endl;
  for (i =9 ; i < n ; i+=9) 
  {
    t = ta + i * h;
    error = exp(t) - 3.0 * cos(t) - y[i];
    outfile1 << setw(15) << t << setw(15) << y[i] << setw(15) << error << endl;
  }
  error = exp(tb) - 3.0 * cos(tb) - beta;
  outfile1 << setw(15) << tb << setw(15) << beta << setw(15) << error << endl;

// BVP2

  const int  m = 4;


  Array x1(m+1,0),x2(n+1,0),x3(n+1,0);
  runtype p, q;

  x1[0] = 1;
  x1[1] = alpha;
  x1[2] = 0;
  x1[3] = alpha;
  x1[4] = 1;
  h = (tb - ta)/n;

  for (i = 1; i <= n; i++)
    {
      rk4sys2(m, h, x1, 1);
      x2[i] = x1[1];
      x3[i] = x1[3];
    }
  p = (beta - x3[n])/(x2[n] - x3[n]);
  q = 1 - p;
  for (i = 1; i <= n; i++)
    x2[i] = p*x2[i] + q*x3[i];
  error = exp(ta) - 3*cos(ta) - alpha;
  outfile2 << setw(15) << "t-Value" << setw(15) << "Solution" << setw(15) << "Error" << endl;
  outfile2 << setw(15) << ta << setw(15) << alpha << setw(15) << error << endl;
  for (i = 9; i < n; i += 9)
    {
      t = ta + i*h;
      error = exp(t) - 3*cos(t) - x2[i];
	  outfile2 << setw(15) << t << setw(15) << x2[i] << setw(15) << error << endl;
    }
  t = ta + i*h;
  error = exp(t) - 3*cos(t) - x2[i];
  outfile2 << setw(15) << t << setw(15) << x2[i] << setw(15) << error << endl;

}