#ifndef __SOLVE4859__
#define __SOLVE4859__
#include "Vector.hh"
#include "Matrix.hh"
#include <assert.h>

enum SolverType {RK4, RK7};

class Solve{
  void (* vfield)(const Vector& x, Vector& y);//vector field
  void (* jacobian)(const Vector& x, Matrix & M);//its Jacobian
  size_t d; //dimension of vector field 

  size_t MaxNumberOfSteps;//maximum number of timesteps allowed
  size_t NumberOfStepsTaken;//number of timesteps taken
  double step; //timestep, can be negative

  int initflag; //initflag=true if soln[0] has been set
  SolverType solver; //solver to be used
  
  Vector *soln; //array of solutions
  Vector solndataspace;//space to store soln
    
  void rk4step(const Vector&, double, Vector&);
  void rk7step(const Vector&, double, Vector&);
  void ex_rk4step(const Vector&, double, Vector&);//ex=extended
  void ex_rk7step(const Vector&, double, Vector&);
  inline void ex_vfield(const Vector&, Vector&);
  
public:
  Solve(void (*vf)(const Vector& , Vector& ),
	void (*jac)(const Vector& , Matrix& ),
	size_t dimension,
	int maxsteps, SolverType s=RK4) ;
  ~Solve();

  void setStart(const Vector& v, double h); //soln[0] = v
  void setStart(const Vector& v){setStart(v,step);};
  void advance(int n=1);//advance soln by n steps

  //methods to get info about solution
  void getSMatrix(double t, Matrix& M); //sensitivity at time t to soln[0]
  void getSMatrix(Matrix& M){ //sensitivity of soln[last step] to soln[0]
    getSMatrix(step*NumberOfStepsTaken,M);
  }
  void shadowSolnDataSpace(Vector& v) const{
    v.shadow(solndataspace);
  }
  void getSoln(int m, Vector& v) const;//v = soln[n]
  void getSoln(int m, int n, Vector& v) const;//v = soln[m..n]
  void getSolnAtT(double t, Vector& v);//v = soln at time t
  size_t getNumberOfStepsTaken() const{
    return NumberOfStepsTaken;
  }
};
       
#endif