Solving Laplace's equation in MATLAB

Question

I've been looking at solving this equation with u(x,0)=u(x,1)=u(0,y)=0 and u(1,y)=sin(pi*y) where (x,y)∈(0,1)x(0,1)

I have solved by hand for h=1/2 and h=1/3 (using the 5-point stencil) however I am struggling with the coding

This is for h=1/3

clear; clc % Clear workspace/command window
% Inputs
Nx = 3;
Ny = 3;
Lx = 1;
Ly = 1;
% Grid
x = linspace (0, Lx, Nx); % Mesh
y = linspace (0, Ly, Ny);
dx = x(2)-x(1); % Mesh size
dy = y(2)-y(1);

% Initialise Matrices

N = 4; % Number of unknowns
A=zeros(N,N); % N rows, N columns
b=zeros(N,1); % N rows, 1 column

% Interior points (not on boundary)

for i=2:Nx-1 % Loop over x direction, first and last grid points skipped
    for j=2:Ny-1 % Loop over y direction, first and last grid points skipped
        n=i+(j-1)*Nx; % Convert i,j grid point to the nth grid point
        A(n,n)=-4; % Main diagonal
        A(n,n-1)=1; % Off diagonal to left
        A(n,n+1)=1; % Off diagonal to right
        A(n,n-Nx)=1; % Far off diagonal to left
        A(n,n+Nx)=1; % Far off diagonal to right
        b(n,1)=0; % Source term
    end
end

% Boundary conditions

% Left boundary condition u = 0

i=1;
for j=1:Ny
    n=i+(j-1)*Nx; % nth row for this ij
    A(n,n)=1; % Main diagonal
    b(n,1)=0; % Boundary condition at y_j
end

% Right boundary condition u = sin(pi*y)

i=N;
for j=1:Ny
    n=i+(j-1)*Nx; % nth row for this ij
    A(n,n)=1; % Main diagonal
    b(n,1)=sin(pi*y(j)); % Boundary condition at y_j
end

% Bottom boundary condition u = 0

j=1;
for i=1:Nx
    n=i+(j-1)*Ny; % nth row for this ij
    A(n,n)=1; % Main diagonal
    b(n,1)=0; % Boundary condition at x_i
end

% Top boundary condition u = 0

j=1;
for i=1:Nx
    n=i+(j-1)*Ny; % nth row for this ij
    A(n,n)=1; % Main diagonal
    b(n,1)=0; % Boundary condition at x_i
end

% Solving A*u=b
% u = inv(A)*b 
u_vec = A\b;

Convert vector solution to a 2D array

 for i=1:Nx
     for j=1:Ny
         n=i+(j-1)*Nx; % nth row for this ij
         u(i,j) = u_vec(n)
     end
 end

 % Plot result
 
surf(x,y,u') % ' takes transpose
xlabel('x')
ylabel('y')

with results of:

u =

NaN

u =

NaN NaN

u =

NaN NaN NaN

u =

NaN NaN NaN NaN 0 0

u =

NaN NaN NaN NaN NaN 0

u =

NaN NaN NaN NaN NaN NaN

u =

NaN NaN NaN NaN NaN NaN NaN 0 0

u =

NaN NaN NaN NaN NaN NaN NaN NaN 0

u =

NaN NaN NaN NaN NaN NaN NaN NaN NaN

any help would be great