Solve Poission equation using Conjugate gradient solver.

I am trying to solve Poisson equation using Conjugate gradient solve.

But the numerical solution does not converge to exact solution....

close all;
N = 100;  % points in the each direction
x0 = 0;
xN = 2;
y0 = 0;
yN = 1;
dx = (xN-x0)/N;
dy = (yN-y0)/N;
dx2_dy2 = dx^2 / dy^2;
x = linspace(x0,xN,N+1);
y = linspace(y0,yN,N+1);
[X Y] = meshgrid(x,y);
Q = 5 * pi^2 * sin(pi .*X) .* sin(2*pi.*Y);
Qint = Q(2:N,2:N);  % obtain interior
q = Qint(:);        % transform into a vector
% True solution
Pd = sin(pi .* X) .* sin(2*pi .*Y);
tol = 1e-5;

Pcg     = zeros(N+1, N+1);  
x       = zeros((N-1)^2, 1); 
x_old   = x;                 
err     = 1;                 
iter    = 0;                 
TN = gallery('tridiag', N-1, -1, 2, -1);
A = kron(speye(N-1), TN) + kron(TN, eye(N-1));
A = A / dx^2;  
b = reshape(Q(2:N, 2:N), [], 1);
r = b - A * x;  
p = r;          
tic;
while (err > tol)
    iter = iter + 1;
    Ap = A * p;
    alpha = (r' * r) / (p' * Ap);
    x_new = x + alpha * p;        
    r_new = r - alpha * Ap;        
    residual_norm = norm(r_new); 
    
    err = norm(x_new - x) 
    if err < tol
        break;
    end
    beta = (r_new' * r_new) / (r' * r);
    p = r_new + beta * p;
    r = r_new;
    x = x_new; 
    
end

err = 125.0411

err = 3.1118e-14

toc;

Elapsed time is 0.074160 seconds.

Pcg(2:N, 2:N) = reshape(x, N-1, N-1);
cg_iter = iter

cg_iter = 2

% Filled Contour Plots
figure(1);
% Numerical Solution Pj
contourf(X, Y, Pcg, 20, 'LineStyle', 'none');  
colorbar;
title('Numerical Solution (Pcg)');
xlabel('x');
ylabel('y');

% error Solution Pj - Pd
figure(2);
contourf(X, Y, Pcg - Pd, 20, 'LineStyle', 'none'); 
title('Error (Pcg)');
xlabel('x');
ylabel('y');
colorbar

[Torsten]

Your matrix A is only suited for an equidistant grid in both directions (x and y).

In your case, it has to be modified because dx = 2*dy.