Я пытаюсь понять, почему размеры этого расчета неверны?

Я решаю уравнение Пуассона и хочу построить график ошибки точного решения в зависимости от количества точек сетки. мой код:

function [Ntot,err] = poisson(N)


nx = N;                % Number of steps in space(x)
ny = N;                % Number of steps in space(y)       
Ntot = nx*ny;
niter = 1000;           % Number of iterations 
dx = 2/(nx-1);          % Width of space step(x)
dy = 2/(ny-1);          % Width of space step(y)
x = -1:dx:1;             % Range of x(-1,1) 
y = -1:dy:1;             % Range of y(-1,1) 
b = zeros(nx,ny);       
dn = zeros(nx,ny);      


% Initial Conditions
d = zeros(nx,ny);                  
u = zeros(nx,ny);

% Boundary conditions
d(:,1) = 0;
d(:,ny) = 0;
d(1,:) = 0;                  
d(nx,:) = 0;


% Source term
b(round(ny/4),round(nx/4)) = 3000;
b(round(ny*3/4),round(nx*3/4)) = -3000;


i = 2:nx-1;
j = 2:ny-1;

% 5-point difference (Explicit)
for it = 1:niter
    dn = d;
    d(i,j) = ((dy^2*(dn(i + 1,j) + dn(i - 1,j))) + (dx^2*(dn(i,j + 1) + dn(i,j - 1))) - (b(i,j)*dx^2*dy*2))/(2*(dx^2 + dy^2));
    u(i,j) = 2*pi*pi*sin(pi*i).*sin(pi*j);

    % Boundary conditions 
    d(:,1) = 0;
    d(:,ny) = 0;
    d(1,:) = 0;                  
    d(nx,:) = 0;
end


%     
% 
% err = abs(u - d);

ошибка, которую я получаю:

Несоответствие размера назначения в подписке.

Ошибка в пуассоне (строка 39)

u(i,j) = 2*pi*pi*sin(pi*i).*sin(pi*j);

Я не уверен, почему он не вычисляет u в каждой точке сетки. Я пытался вывести его из цикла for, но это не помогло. Любые идеи были бы хорошы.