/*
*    Euler's method  to 
*    solve y' = y-t^2+1, y(0) = 0.5, 0 <= t <= 2.
*
*    The exact solution is y(t) = (t+1)^2 - 0.5e^t   
*    INPUT:   ENDPOINTS A,B; INITIAL CONDITION ALPHA; INTEGER N.
*
*    OUTPUT:  APPROXIMATION W TO Y AT THE (N+1) VALUES OF T.
*/

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

#define true 1
#define false 0
#define MAX_N 2000

static   double F(double, double);
static   double exact_soln(double); 

using namespace std;

int main()
{
      double A,B,ALPHA,H,T,W[MAX_N], K1,K2,K3,K4;
      int I,N;

      A = 0.0;
      B = 2.0;
      N = 10;

      cout.setf(ios::fixed,ios::floatfield);
      cout.precision(10);

      /* STEP 1 */
      H = (B - A) / N;
      T = A;
      W[0] = 0.5; // initial value

      /* STEP 2 --- Use Euler's method */
      /// NOTE: The "for" loop starts with I = 1. 
      for (I=1; I<=N; I++) 
      {
         /* COMPUTE W(I) */
         W[I] = W[I-1] + H*F(T, W[I-1]);

         /* COMPUTE T(I) */ 
         T = A + I * H;

         /* STEP 5 */
         cout <<"At time "<< T <<" numerical solution = "<< W[I] <<"; exact soln = "<< exact_soln(T) << endl;
      }

      return 0;
}
/*  Change function F for a new problem   */
double F(double T, double Y)
{
   double f; 

   f = Y - T*T + 1.0;
   return f;
}

double exact_soln(double T)
{
    return ((T+1.0)*(T+1.0)- 0.5*exp(T));
}