% Full multigrid V-cycle on u'' = f(x) = 6x, u(0) = a = 0, u(1) = b = 1 for 
% h = 1/32

  a = 0;
  b = 1;

  nu1 = 3;                 % Pre- and post-relaxation iterations
  nu2 = 1;

  pmax = 3;                % Maximum level of refinement h = 2^{-pmax}
  N = 2^pmax;
  h = 1/N;
  x = h*(0:N);
  f = 6*x;

  r = zeros(1,N+1);
  u = zeros(1,N+1);
  u(1) = a;
  u(N+1) = b;
  r(2:N) = f(2:N) - (u(1:N-1) - 2*u(2:N) + u(3:N+1))/h^2;  % residual
  normf = norm(r);      % 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.

%  Plot exact solution

  plot(x,x.^3,'r')
  xlabel('x')
  ylabel('u')
  title(['FMV, h = 1/8, \nu_1=' num2str(nu1) ', \nu_2=' num2str(nu2)])
  hold on
  char = '*x+o.';

%  Full multigrid solve


%  Exact solve on coarsest grid h = 1/2

  hp = 1/2;
  u(1) = a;
  u(3) = b;
  u(2) = 0.5*(u(1) + u(3) - hp^2*f(N/2+1));
  plot(0:0.5:1,u(1:3),char(1))
  text(0.1,0.85,[char(1) ' h =1/2'])

  for p = 2:pmax
    Np = 2^p;
    hp = 1/Np;

%    Interpolate to get u from next-coarser grid

    v = u(1:(Np/2+1));
    u = zeros(1,Np+1);
    u(1) = v(1);
    j = 1:(Np/2);
    twoj = 2*j;
    twojp1 = 2*j+1;
    u(twoj) = (v(j) + v(j+1))/2;
    u(twojp1)= v(j+1);

%    Restrict (via injection) from finest-grid f to get fp at level p

    stride = N/Np;
    jp = 1:stride:(N+1);
    fp = f(jp);

%  Do a V-cycle.

    u = vcycle(u,fp,hp,nu1,nu2);
    plot(x(jp),u,char(p))
    text(0.1,1-0.15*p,[char(p) ' h = 1/' num2str(Np)])
  end
  hold off

%  Calculate residual vector and relative residual after FMV

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


