#include <iostream>
#include <math.h>
#include "pstools.h"

using namespace std;

int main()
{
  int order = 350;
  int numEqs = 3;
  double c[numEqs][order+1];
  double endPt = 1.0;
  double stepSize = 0.01;
  int numSteps = endPt / stepSize;
  
  double v = 10.4; //must be more than 1
  double k = 5.5;
  
  c[0][0] = 0.0;
  c[1][0] = 1.0/v;
  c[2][0] = 1.0;
  

  double trajectory[numSteps];
  trajectory[0] = c[0][0];
  c[2][1] = 1;//?
  
  for(int step=0;step<=numSteps;step++)
  {
    
    for(int i=0;i<=order;i++)
    {
       c[0][i+1] = cp(c[1],c[2],i)/(i+1);
       c[1][i+1] = -cp(c[1],c[1],i)/(i+1);
       if(i>0)
         c[2][i+1] = -2.0*k*c[2][i-1]/(i+1);
    }
    
    #pragma omp parallel for
    for(int i=0;i<3;i++)
    {
      c[i][0]=polyval(c[i],stepSize,order);
    }
    trajectory[step+1]=c[0][0];
  }
  
  //cout<<c[0][0]<<"\n";
  for(int i = 0;i<=5;i++)
  {
    cout<<trajectory[i]<<" ";
  }
  cout<<endl;

  
  return 0;
}