double f(double x,double v,double t)
/*
    x position
    v velocity
    t time
*/
{
  double k=1.0; // spring.
  return -k*x;
}

void IntegrateEL(double dt, double& x, double& v, double& t)
// Euler method  f is a function that evalues the force for a given x,v,t
{
    double a=f(x,v,t)/m;

    double xnext=x+v*dt;
    double vnext=v+a*dt;

    x=xnext;
    v=vnext;
    t+=dt;
}


void IntegrateMP(double dt,  double& x, double v&, double& t)
// Midpoint method  f is a function that evalues the force for a given x,v,t
{
    double xm=x+0.5*dt*v;
    double vm=v+f(x,v,t)*0.5*dt/m;

    double xnext = x + v*dt + f(x,v,t)*0.5*dt*dt/m;
    double vnext = v + f(xm, vm, t+0.5*dt)*dt/m;

    x=xnext;
    v=vnext;
    t+=dt;
}

void IntegrateRK4(double dt,  double& x, double& v, double& t)
// Runge-Kutta RK4
{
// derivatives at t
   double f1x=v;
   double f1v=f(x,v,t)/m;

// First estimate of derivatives at t+=dt/2
   double f2x=v+ 0.5*dt*f1x;
   double f2v=f(x+0.5*dt*f1x,v+0.5*dt*f1v,t+0.5*dt)/m;

// Corrected estimates at t+=dt/2
   double f3x=v+0.5*dt*f2x;
   double f3v=f(x+0.5*dt*f2x,v+0.5*dt*f2v,t+0.5*dt)/m;

// Final estimates at t+=dt
   double f4x=v+dt*f3x;
   double f4v=f(x+dt*f3x,v+dt*f3v,t+dt)/m;

   double xnext=x+dt*(f1x+2.0*f2x+2.0*f3x+f4x)/6.0;
   double vnext=v+dt*(f1v+2.0*f2v+2.0*f3v+f4v)/6.0;

   x=xnext;
   v=vnext;
   t+=dt;
}

void IntegrateVerlet(double dt,  double& x, double& v, double& t)
// Verlet (Velocity version)
{
   double a=f(x,v,t)/m;
   double vhalf=v+0.5*dt*a;
   double xnext=x+dt*vhalf;
   double anext=f(xnext,vhalf,t+dt)/m;
   double vnext=vhalf+0.5*dt*anext;

   x=xnext;
   v=vnext;
   t+=dt;
}