% Multigrid V-cycles on u'' = f(x) = 6x, u(0) = a = 0, u(1) = b = 1 for h = 1/8

  a = 0;
  b = 1;

  maxiter = 10;             % Number of V-cycles
  nu1 = 1;                 % Number of pre-relaxations at each level.
  nu2 = 1;                 % Number of post-relaxations at each level.

  pmax = 3;                % Maximum level of refinement h = 2^{-pmax}
  N = 2^pmax;
  h = 1/N;
  x = h*(0:N);
  f = 6*x;
  u = [a zeros(1,N-1) b];  % Initial guess is zero at all interior gridpoints

  r = zeros(1,N+1);
  r(2:N) = f(2:N) - (u(1:N-1) - 2*u(2:N) + u(3:N+1))/h^2;  % Initial residual
  normf = norm(r(2:N))  % If problem is phrased as Au=f, norm(f) is equal
                        % to the norm of the residual computed above, since
                        % initially u = 0 at all interior gridpoints.

%  Do maxiter V-cycles.

  subplot(2,2,1)
  plot(x,u,'gx')
  xlabel('x')
  ylabel('u')
  title('Multigrid V-Cycles')
  hold on
  plot(x,x.^3,'r')
  for k = 1:maxiter
    u = vcycle(u,f,h,nu1,nu2);
    if (k < 4)
      plot(x,u,'b+')
    end

%  Calculate residual vector and relative residual after each V-cycle


    r(2:N) = f(2:N) - (u(1:N-1) - 2*u(2:N) + u(3:N+1))/h^2;  % residual
    relres(k) = norm(r)/normf;  % relative residual
  end
  hold off

  subplot(2,2,2)
  semilogy(0:maxiter, [1 relres(1:maxiter)],'+')
  axis([0 maxiter 1e-6 1])
  xlabel('Iteration')
  ylabel('Relative Residual')
  title(['h = 1/' num2str(N)])

