#include <iostream>
#include <string>
#include <cmath>
using namespace std;

#define N 20

double fex (double t);
double f0 (double t, double y);
double f1 (double t, double y, double y1);
double f2 (double t, double y, double y1, double y2);
double f3 (double t, double y, double y1, double y2, double y3);

int main()
{
  double T = 1.0;
  double y0 = 1.0;
  double h, h2, h3, h4;
  double H[4], t[N], yex[N], y[N], c[N][4], err[N];
  double tt, yy0, yy1, yy2, yy3, yy4, sum;

  h = T/N;
  h2 = h*h/2.0;
  h3 = h2*h/3.0;
  h4 = h2*h2/6.0;
  H[0]=1.0; H[1]=h; H[2]=h2; H[3]=h3; H[4]=h4;

  t[0] = 0.0;
  yex[0] = y[0] = y0;
  c[0][0] = y0;
  c[0][1] = f0(t[0],y[0]);
  c[0][2] = f1(t[0],y[0],c[0][1]);
  c[0][3] = f2(t[0],y[0],c[0][1],c[0][2]);
  c[0][4] = f3(t[0],y[0],c[0][1],c[0][2],c[0][3]);
  err[0] = 0.0;
 
  cout << "t[n]\t    y[n]\t c[n-1][1]    c[n-1][2]\t   yex[n]\t err[n]" << endl;
  for (int j=0;j<72;j++) cout << "-"; 
  cout << endl;

  printf("%5.3f    %12.8f  %8.4f   %8.4f", t[0], y[0],  0.000,  -1.000);
  printf("\t %12.8f    %10.6e\n", yex[0], err[0]);

  for(int j=0;j<N;j++) {
    tt = t[j];
    yy0 = y[j];
    c[j][0] = yy0;
    yy1 = f0(tt, yy0); 
    c[j][1] = yy1;
    yy2 = f1(tt,yy0,yy1);
    c[j][2] = yy2;
    yy3 = f2(tt,yy0,yy1,yy2);
    c[j][3] = yy3;
    yy4 = f3(tt,yy0,yy1,yy2,yy3);
    c[j][4] = yy4;
           // the loop lets itself do more effectively
    sum = 0.0;
    for(int k=0;k<5;k++)
        sum = sum + c[j][k]*H[k];

    y[j+1] = sum;
    t[j+1] = tt + h;
    yex[j+1] = fex(tt+h);
    err[j+1] = y[j+1] - yex[j+1];

    printf("%5.3f    %12.8f  %8.4f   %8.4f",
                             t[j+1], y[j+1],  c[j][1], c[j][2]);
    printf("\t %12.8f    %10.6e\n", yex[j+1], err[j+1]);
  }
}
  
// ======================================================

double f0 (double t, double y)
{  
    // RHS of f(t,y(t))
    return (-2*t*y*y);
}

double f1 (double t, double y, double y1)
{
    // RHS of d/dt(f(t,y(t))) 
    return (-2*y*y-4*t*y*y1);
}

double f2 (double t, double y, double y1, double y2)
{
    // RHS of d^2/dt^2(f(t,y(t))) 
    return ( -8*y*y1  -4*t*( y1*y1 + y*y2 ));
}

double f3 (double t, double y, double y1, double y2, double y3)
{
    // RHS of d^3/dt^3(f(t,y(t)))
    return ( -12*(y1*y1 + t*y1*y2 + y*y2) - 4*t*y*y3 ); 
} 

// ======================================================

double fex (double t)
{
   // analytic solution
   return ( 1.0/(1.0 + t*t) );
} 

// ======================================================