////////////////////////////////////////////////////////////////////////
void TridiagonalSolve(double *am, double *a0, double *ap, double *b,
                      double *phi, int n)
//
// Use recursion formula to solve a tridiagonal system of the type
//
//  am[i] \phi[i-1] + a0[i] \phi[i] + ap[i] \phi[i+1] = b[i]
//
// from i = 1..n-1.
// Boundary coefficients must have been imposed at \phi[0] and \phi[n]
// On output \phi[] is replaced with the updated solution.
////////////////////////////////////////////////////////////////////////
{
  using namespace std;
  int i;
  double alpha[n], beta[n], gamma;


  alpha[n-1] = 0.0;
  beta[n-1]  = phi[n];
  for (i = n-1; i > 0; i--){
    gamma      = -1.0/(a0[i] +   ap[i]*alpha[i]);
    alpha[i-1] = gamma*am[i];
    beta[i-1]  = gamma*(ap[i]*beta[i] - b[i]);
  }

  for (i = 0; i < n-1; i++){
    phi[i+1] = alpha[i]*phi[i] + beta[i];
  }

// Verify solution of tridiagonal matrix

  double res;
  for (i = 1; i < n; i++){
    res = am[i]*phi[i-1] + a0[i]*phi[i] + ap[i]*phi[i+1] - b[i];
    if (fabs(res) > 1.e-9){
      cout << "! TridiagonalSolve: not correct" << endl;
      exit(1);
    }
  }
}