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
Afun1=@(v) Mfun1(Ahat*v,A,nu);   % action of the preconditioned matrix 1
Mfhat1=Mfun1(fhat,A,nu);         % rhs of the preconditioned system 1
u0=rand(3*m^d,1);                % random initial guess
kmax=min(200,m^d);               % set max number of iterations
[u1,~,~,~,rv1]=gmres(Afun1,Mfhat1,[],1e-6,kmax,[],[],u0);
At=spdiags(diag(A),0,m^d,m^d);
Afun2=@(v) Mfun2(Ahat*v,At,nu);  % action of the preconditioned matrix 2
Mfhat2=Mfun2(fhat,At,nu);        % rhs of the preconditioned system 2
[u2,~,~,~,rv2]=gmres(Afun2,Mfhat2,[],1e-6,kmax,[],[],u0);
semilogy(0:length(rv1)-1,rv1./rv1(1),'b'); hold on;
semilogy(0:length(rv2)-1,rv2./rv2(1),'r');
axis([0 40 1e-6 1]); yticks([1e-6 1e-4 1e-2 1]);
legend('Schur1','Schur2','location','NorthEast');
xlabel('iterations');
ylabel('norm of the residual'); 
set(gca,'fontsize',18,'linewidth',2);