For dx = dy = 0.01, Nx = Ny = 101, not 30 in your code. I just realized this after setting up the code below.
dx = 0.01;
dy = 0.01;
x = 0:dx:1;
y = 0:dy:1;
nx = numel(x);
ny = numel(y);
a =1/dx^2;
c =1/dy^2;
b =-2*(a+c);
A = zeros(nx*ny,nx*ny);
B = zeros(nx*ny,1);
% Boundaries
% Boundary values at y = 0
for ix = 1:nx
A(ix,ix) = 1.0;
B(ix) = 0.0;
end
% Boundary values at x = 0
for iy = 2:ny-1
k = nx*(iy-1) + 1;
A(k,k) = 1.0;
B(k) = 0.0;
end
% Boundary values at x = 1
for iy = 2:ny-1
k = nx*iy;
A(k,k) = 1.0;
B(k) = y(iy);
end
% Boundary values at y = 1
for ix = 1:nx
k = nx*(ny-1) + ix;
A(k,k) = 1.0;
B(k) = x(ix);
end
% Inner grid points
for iy = 2:ny-1
for ix = 2:nx-1
k = (iy-1)*nx + ix;
A(k,k) = b;
A(k,k+1) = a;
A(k,k-1) = a;
A(k,k+nx) = c;
A(k,k-nx) = c;
end
end
u = A\B;
%u
U = zeros(nx,ny);
for iy = 1:ny
for ix = 1:nx
k = (iy-1)*nx + ix;
U(ix,iy) = u(k);
end
end
[X,Y] = meshgrid(x,y);
surf(X, Y, U);
