function u = vcycle(u,f,h,nu1,nu2)

%  Does multigrid V-cycle on FDA for u'' = f, with nu1 pre-relaxation and
%  nu2 post-relaxation sweeps at 
%  each level, and returns an updated u.  Note that the first and last values
%  of u and f are the values at the boundaries x = 0 and 1.

  om = 2/3;
  N = length(u) - 1;
  if (N == 2)

%  Calculate exact soln.

    u(2) = 0.5*(u(1) + u(3) - h^2*f(2));
    return
  else

%  weighted Jacobi pre-relaxation

    int = 2:N;    % Interior indices 
    for nu = 1:nu1
      u(int) = om*0.5*(u(int-1) + u(int+1) - h^2*f(int)) + (1 - om)*u(int);
    end

%  Calculate residual on coarse grid points (i. e. restrict r via injection)

    r = zeros(1,N/2+1);
    int2 = 3:2:(N-1);  % Interior grid points that are also on 2h grid
    r(2:(N/2)) = f(int2) - (u(int2-1) - 2*u(int2) + u(int2+1))/h^2;

%    Initialize and calculate coarse grid correction

    v = zeros(1,N/2+1);
    v = vcycle(v,r,2*h,nu1,nu2);

%    Interpolate coarse grid correction and add to u  

    intcoarse = 2:(N/2);     % Interior 2h-grid point index range.
    intnot =[int2-1 N];      % Interior h-grid points not on 2h grid
    u(int2)= u(int2) + v(intcoarse);
    u(intnot) = u(intnot) + 0.5*(v(1:(N/2)) + v(2:(N/2+1)));

% Post-relaxation on corrected solution

    for nu = 1:nu2
      u(int) = om*0.5*(u(int-1) + u(int+1) - h^2*f(int)) + (1 - om)*u(int);
    end
  end
