#include <stdio.h>   //standard IO
#include <string.h>  //string manipulations
#include <math.h>    // alternate math functions

//Define a global parameter
const double accel = 4.0 ;

//Dependences of derivative function
//x is independent variable
//y is dependent vector
//dydx is dy/dx vector at x, y
void Dydx(double x, double y[], double dydx[]) ;

//Dependences of step function
//n is the number of differential equations
//dx is stepsize
//*derivs is the derivative function (Dydx above)
void Step(double x, double y[], double dydx[], int n, double dx,
           void (*derivs)(double, double [], double [])) ;

//main program
int main()
{
const int nphys = 2 ;
int ndim ;
double y[nphys], dydx[nphys], x, t, dt, t0, tf, x0, v0, yexact[nphys] ;

FILE *output ; output = fopen("y_vs_x.dat","w") ;
FILE *outexa ; outexa = fopen("y_vs_x_exact.dat","w") ;

//initialize t0, tf = t_final, dt, x0 & v0
t0 = 2.1 ; tf = 5.7 ; dt = 0.01 ; x0 = 16.0 ; v0 = -7.3 ;

//fill the array containing the y_vector at the initial time
//ndim is the number of elements in y
y[0] = x0 ; y[1] = v0 ; t = t0 ; ndim = 2 ;

//the following while loop steps through time and outputs the results
//to two files; one file contains the exact results, the other the
//numerical approximation
  while(t < tf)
  {
//step t then step y
  t += dt ;
  Step(t,y,dydx,ndim,dt,Dydx) ;

//compute the exact value of y for constant acceleration
  yexact[0] = x0 + v0*(t-t0) + 0.5*accel*(t-t0)*(t-t0) ;
  yexact[1] = v0 + accel*(t-t0) ;

//print out the numerical approximation and the exact value
  fprintf(output,"%13.6E %13.6E %13.6E\n",t,y[0],y[1]) ;
  fprintf(outexa,"%13.6E %13.6E %13.6E\n",t,yexact[0],yexact[1]) ;
  }
}


//Function for the calculating the derivative
void Dydx(double x, double y[],double dydx[])
{
int j ;
//dx/dt = v
dydx[0] = y[1] ;
//d^2 x/dt^2 = Force
dydx[1] = accel ;
}


//Step moves y forward by dx using the lowest order approximation
//(known as the Euler method)
void Step(double x, double y[], double dydx[], int n, double dx,
           void (*derivs)(double, double [], double []))
{
 int j ;
//call the function that computes the derivative of y
 (*derivs)(x,y,dydx) ;
 for(j = 0 ; j < n ; j++) y[j] += dydx[j]*dx ;
}