m=15;                          % number of mesh points in each direction
d=2;                           % dimension of the problem
nu=1e-3;                       % regularization parameter
h=1/(m+1);                     % mesh size
A=Laplacian(m,d)/h^2;          % negative 2D discrete Laplacian
Ahat=[sparse(m^d,m^d), A, -speye(m^d); ...
      A, speye(m^d), sparse(m^d,m^d); ...
      -speye(m^d), sparse(m^d,m^d), nu*speye(m^d)]; % all-at-once operator
f=-ones(m^d,1);                % discrete f=-1
yd=ones(m^d,1);                % discrete y_d=1
fhat=[f;yd;zeros(m^d,1)];      % assemble the right-hand side vector
u0=rand(3*m^d,1);              % random initial guess
[u,res] = MG_OC(Ahat,fhat,nu,u0,3,1e-16,21);    
semilogy(0:20,res(1:21)./res(1),'o-b');
axis([0 20 1e-10 1]); yticks([1e-10 1e-8 1e-6 1e-4 1e-2 1]);
xlabel('iterations');
ylabel('norm of the residual'); 
set(gca,'fontsize',18,'linewidth',2);
