/*
*     TRAPEZOIDAL WITH NEWTON ITERATION ALGORITHM 5.8
*
*     TO APPROXIMATE THE SOLUTION OF THE INITIAL VALUE PROBLEM:
*                Y' = F(T,Y), A <= T <= B, Y(A) = ALPHA,
*     AT (N+1) EQUALLY SPACED NUMBERS IN THE INTERVAL [A,B].
*
*     INPUT:   ENDPOINTS A,B; INITIAL CONDITION ALPHA; INTEGER N;
*              TOLERANCE TOL; MAXIMUM NUMBER OF ITERATIONS M AT ANY ONE STEP.
*
*     OUTPUT:  APPROXIMATION W TO Y AT THE (N+1) VALUES OF T
*              OR A MESSAGE OF FAILURE.
*/

#include<stdio.h>
#include<math.h>
#include<iostream> 
#define true 1
#define false 0

static   double F(double, double);
static   double FYP(double, double);
static   double exact_soln(double T);
static   double absval(double);

using namespace std; 

int main()
{
      double A,B,ALPHA,TOL,W,T,H,XK1,W0,Y;
      int N,M,I,J,IFLAG,OK;

      A = 0.0;
      B = 1.0;
      ALPHA = -1.0;
      N = 10;

      /* For Newton's method */
      TOL = 1.0e-12;  // tolerance for testing convergence 
      M   = 100;      // max. number of iteration

      /* STEP 1 */
      W = ALPHA;
      T = A;
      H = (B - A) / N;
      I = 1;
      OK = true;

      /* STEP 2 */
      for (I = 1; I <= N; I++) 
      {
         /* STEP 3, Newton's method to find W_I */ 
         XK1 = W + 0.5 * H * F(T, W);
         W0 = XK1;
         J = 1;
         IFLAG = 0;
         /* STEP 4 */
         while ((IFLAG == 0) && OK) 
         {
            /* STEP 5, Newton's method */
            W = W0 - (W0 - XK1 - 0.5 * H * F(T + H, W0)) /
                (1.0 - 0.5 * H * FYP( T + H, W0));
            /* STEP 6 */
            if (absval(W-W0) < TOL) 
            {
               IFLAG = 1;
               /* STEP 7 */
               T = A + I * H;
               cout << "At time " << T <<", Num. soln = " << W <<" Exact soln " <<exact_soln(T)  << ". Num of iter. " << J <<endl;  
            }  
            else 
            {
               J++;
               W0 = W;
               if (J > M) OK = false;    
            }   
         }  
         
         /* Netwon's method failed, terminate the calculation */
         if (!OK)
         {
            cout <<"Maximum Number of Iterations Exceeded, Newton's method failed" << endl;
            break;
         }  
      }
      /* STEP 8 */

   return 0;
}

/*  Change functions F and FYP for a new problem.   */
static double F(double T, double Y)
{
   double f; 

   f = 5*exp(5*T)*(Y-T)*(Y-T)+1;
   return f;
}
/*  Function FYP is the partial derivative of F with respect to y.  */
static double FYP(double T, double Y)                         
{
   double f; 

   f = 10*exp(5*T)*(Y-T); 
   return f;
}


/* Absolute Value Function */
static double absval(double val)
{
   if (val >= 0) return val;
   else return -val;
}

static double exact_soln(double T)
{
    return (T - exp(-5.0*T));
}