// FILE DemoStiffDetestD1.cc

#include "VNODE.h"
#include "StiffDetestD1.h"


int main()
{
  VnodeInit();
  
  PtrODENumeric ODE = new StiffDetestD1; 
  PtrDataRepr DataRepr = new TAYLOR_EXPANSION<StiffDetestD1TaylGenODE,StiffDetestD1TaylGenVAR>;
  
  ODE->LoadProblemParam(1);
  
  PtrODESolver Solver = new SOLVER_2( ODE, DataRepr );
  
  cout << "*** Integrating " << ODE->GetName() << endl;

  double start = 100;
  double end = 400;
  double step = 100;

  for ( double t = start; t <= end; t += step )
    {
      // If t != start, we continue the integration from the current point.
      Solver -> IntegrateTo( t, t != start );
      
      cout << "\n   Enclosure at t = " << ODE->GetTcur() << endl;
      PrintVec( ODE->GetTightEncl() );
    }

  // We save the solution at the end point.
  INTERVAL_VECTOR Yend( ODE->GetTightEncl() );
  
  // Integrate the same problem from ODE->GetT0() to the endpoint.
  Solver->IntegrateTo( end, false );
  cout << "\n   Enclosure at t = " << ODE->GetTcur() << endl;
  PrintVec( ODE->GetTightEncl() );

  /* Yend and ODE->GetTightEncl() must intersect. However, they may
     not be the same for the following reasons.  

     The stepsize controller may choose a much smaller stepsize, when
     the integration is close to an end point. In the for loop, we
     have more then one end points, while in the integration outside
     the for loop, there is only one end point. Thus, the number of
     steps may be different, and the enclosures at the point "end" may
     be slightly different.  */

  assert( Intersection(Yend, Yend, ODE->GetTightEncl()) );

  return 0;
}