#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 ;

//main program
int main()
{
double x, v, a, t, dt, t0, tf, x0, v0, xexact, vexact ;

FILE *output ; output = fopen("y_vs_x_simp.dat","w") ;
FILE *outexa ; outexa = fopen("y_vs_x_simp_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 ;

//set the initial position, velocity and time
//ndim is the number of elements in y
x = x0 ; v = v0 ; t = t0 ;

//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
  t += dt ;
//compute the acceleration
  a = accel ;
//step x
  x += v*dt ;
//step v
  v += a*dt

//compute the exact value of y for constant acceleration
  xexact = x0 + v0*(t-t0) + 0.5*accel*(t-t0)*(t-t0) ;
  vexact = v0 + accel*(t-t0) ;

//print out the numerical approximation and the exact value
  fprintf(output,"%13.6E %13.6E %13.6E\n",t,x,v) ;
  fprintf(outexa,"%13.6E %13.6E %13.6E\n",t,xexact,vexact) ;
  }
 fclose(output) ;
 fclose(outexa) ;
}