/*
	test harness for ODE integrators
*/

#include <odeint.h>
#include <function.h>

class integrator: public ODEvector, public Classio
{
   public:
   Array<Realfunction1d *> fn;
   
   void Read(istream& s) { s >> fn; }
   void Write(ostream& s) const { s << fn << endl; }
   
   Vector dfdx(const double x, const Vector& f) const
   {
      assert(f.n == fn.n);
      Vector df(f.n); 
      for (int i=0;i<fn.n;i++) df[i] = (*(fn[i]))(x); 
      return df; 
   }
};

int main()
{
integrator A; cin >> A; A.set_defaults();
double x1, x2; cin >> x1 >> x2;
Vector f1(A.fn.n); f1.Read_elements(cin);

int nsteps, hitmax, underflow;
Vector f2 = A.rkaccint(x1,x2,f1,nsteps,hitmax,underflow);
cout << "f(x2) = "; f2.Write_elements(cout); cout << endl;
cout << nsteps << " steps\n";
if (hitmax) cout << "...max steps\n";
if (underflow) cout << "stepsize underflow\n";

return 0;
}