//  Solves linear equations for simple tridiagonal matrix using the iterative Jacobi method
....
// Begin main program
int main(int argc, char *argv[]){
     // missing statements, see code link above

      mat A = zeros<mat>(n,n);
      // Set up arrays for the simple case
      vec b(n);  vec x(n);
      A(0,1) = -1;  x(0) = h;  b(0) =  hh*f(x(0)); 
      x(n-1) = x(0)+(n-1)*h; b(n-1) = hh*f(x(n-1)); 
      for (int i = 1; i < n-1; i++){ 
        x(i) = x(i-1)+h; 
	b(i) = hh*f(x(i));
        A(i,i-1)  = -1.0;
        A(i,i+1)  = -1.0;
      }
     A(n-2,n-1) = -1.0; A(n-1,n-2) = -1.0;
  // solve Ax = b by iteration with a random starting vector
     int maxiter = 100; double diff = 1.0; 
     double epsilon = 1.0e-10;  int iter = 0;
      vec SolutionOld  = randu<vec>(n);
      vec SolutionNew  = zeros<vec>(n);
      while (iter <= maxiter || diff > epsilon){
	SolutionNew = (b -A*SolutionOld)*0.5; 
        iter++; diff = fabs(sum(SolutionNew-SolutionOld)/n);
        SolutionOld = SolutionNew;
      }
      vec solution = SolutionOld;}