N = 2048;
h = 1/N;
x = (0:h:1)';      % N+1 points
f = - 1 + 100*cos(24*pi*x) + 5e4*cos(1e3*pi*x);

% Direct solver
e = ones(N+1,1);
A = spdiags([e -2*e e], -1:1, N-1, N-1)/h^2;
u = A \ f(2:(end-1));
u_exact = [0; u; 0];

% Multigrid
tic,
u = solvepoissonmultigrid(f);
toc
figure(1); clf;
plot(x,u_exact,'b-'); hold on;
plot(x,u,'rx');
disp(['Error = ',num2str(max(abs(u-u_exact)))]);
