void crank_nicolson(int n, int tsteps, double delta_x, double alpha)
{
   double a, b, c;
   vec u(n+1); // This is u in Au = r
   vec r(n+1); // Right side of matrix equation Au=r
   ....
   // setting up the matrix 
   a = c = - alpha;
   b = 2 + 2*alpha;

   // Time iteration
   for (int t = 1; t <= tsteps; t++) {
      // Calculate r for use in tridag, right hand side of the Crank Nicolson method
      for (int i = 1; i < n; i++) {
	 r(i) = alpha*u(i-1) + (2 - 2*alpha)*u(i) + alpha*u(i+1);
      }
      r(0) = 0;
      r(n) = 0;
      //  Then solve the tridiagonal matrix
      tridiag(a, b, c, r, u, xsteps+1);
      u(0) = 0;
      u(n) = 0;
      //  Eventual print statements etc
      ....
}