C++でファン・デル・ポル方程式を後退オイラー法で解く

ファン・デア・ポル方程式を後退オイラー法で解くプログラム

#include <iostream>
#include <math.h>
using namespace std;
 
void back_euler(double&,double&,double,double);
 
int main(){
  double x, v, f, h, t, tend;
  int step;
 
  x = 1.0;
  v = 0.0;
  f = 0;
  tend = 100;
  h = 0.1;
 
  step = (int)(tend/h);
 
  cout << t << " " << x << " " << v << endl;
  for (int i = 0; i < step; i++){
    back_euler(x,v,f,h);
    t += h;
    cout << t << " " << x << " " << v << endl;
  }
}
 
void back_euler(double &x, double &v, double f, double h){
  double dfxdx, dfxdv, dfvdx, dfvdv;
  double idet;
  double fx, fv;
  double xold, vold, delta, e;
  double x0, v0;
 
  x0 = x;
  v0 = v;
 
  do{
  dfxdx = 1.0;
  dfxdv = -h;
  dfvdx = h*(2*f*v*x + 1);
  dfvdv = 1 + h*f*(x*x -1);
 
  fx = x - h*v - x0;
  fv = v + h*(f*v*(x*x -1) + x) - v0;
 
  idet = 1.0/(dfxdx*dfvdv - dfxdv*dfvdx);
 
  xold = x;
  vold = v;
 
  x = x - idet*(dfvdv*fx - dfxdv*fv);
  v = v - idet*(-dfvdx*fx + dfxdx*fv);
 
  delta = sqrt(fx*fx + fv*fv);
  e = sqrt((x-xold)*(x-xold) + (v-vold)*(v-vold));
  }while(delta > 0.001 && e > 0.001);
}