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
Afun=@(v) nu*v+(A\(A\v));        % define function Afun
f=-ones(m^d,1);                  % discrete f=-1
yd=ones(m^d,1);                  % discrete y_d=1
rhs=-A\(A\f-yd);                 % assemble the right-hand side vector
u0=rand(m^d,1);                  % random initial guess
kmax=min(100,m^d);               % set max number of iterations
[u1,~,~,~,rv1]=gmres(Afun,rhs,[],1e-6,kmax,[],[],u0);
[u2,~,~,~,rv2]=pcg(Afun,rhs,1e-6,kmax,[],[],u0);
semilogy(0:length(rv1)-1,rv1./rv1(1),'o-b'); hold on;
semilogy(0:length(rv2)-1,rv2./rv2(1),'x-r');
axis([0 50 1e-8 1]); yticks([1e-8 1e-6 1e-4 1e-2 1]);
legend('GMRES','CG','Location','NorthEast');
xlabel('iterations');
ylabel('norms of the residuals'); 
set(gca,'fontsize',18,'linewidth',2);