#include <stdio.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include "dumka3.cpp"
//#include "euler.cpp"
#include <vector>
//CC HeatEquationNoVector.cpp -lm -LANG:std -64 -apo
using namespace std;
ofstream  solution( "solution" ); // open output file for solution


template <class Vector, class Element>
class VectorFunction: public VectorFunctionBase<Vector,Element> {
protected:
  double m_mu;
public:
  VectorFunction() {
  }
  VectorFunction ( double mu ) {
    this->m_mu = mu;
  }
  bool isStepAccepted;
  void compute(Vector& y, double t, Vector& result) {
    long size = y.size();
    long size1= size - 1;  
    double h=1./(double)size1;
    double h_inv=1./h;
    double h_inv2=m_mu*h_inv*h_inv;
    double h_inv05=0.25*h_inv;
//cout << h << " " << h_inv << " " << h_inv2 << " " << h_inv05 << endl;
    y[0]=0.;
    y[size1]=0.;
    for ( int i = 1; i < size1; i++ ) {  
      double right = y[i+1];
      double left = y[i-1];
      double center = y[i];
      result[i] = -h_inv05*(right*right-left*left) + h_inv2*(right-2.*center+left);
    }
    result[0] = 0.;
    result[size1]=0.;  
    if ( isStepAccepted ) {
      isStepAccepted = false;
      //output(solution, y, t);
    }
  };
  Element GetNonZeroElement() {     
     return (rand()%100)/100.;
   }
  double cour(Vector& y, double t) {
    return 3.53695e-05;//0.019;//0.038948
  };
};
template<class Vector>
void output(ofstream& solution, Vector& y, double t) {
   long size = y.size();
   for ( int i = 0; i < size; i++ ) {
    double x = ((double)i)/((double)size-1.);
    solution << t << " " << x << " " << y[i] << endl; 
  }
}
int main ( int argc, char** argv ) {
  // y - solution of ODE 
  int size = 5001;
  vector<double>  y(size);
  for ( int i = 0; i < size; i++ ) {   
    double x = ((double)i)/((double)size-1.);
    y[i] = 1.5*x*(1.-x)*(1.-x);
  } 
  VectorFunction<vector<double>,double> f( 0.0005 );

  double tinitial=0.;
  double t=tinitial;    // t - initial time
  double tend=2.5;       // tend - interval of integration
  double rtol = 0.1; // rtol - relative error 
  double atol=0.1;      // tolerance of computation
  double h0 = 0.001;  // h0 - initial step size
                      // f - class VectorFunction which has 
  output(solution, y, t);
  dumka3(t, tend, h0, atol, rtol, f, y, true );
  //euler( t, tend, f, y );                    //     method Vector  
  output(solution, y, t);
  solution.close();

  return 1;
}