How can I get more accurate electric scalar potential in 2D closed box?

Question

I am trying to use poisson equation to plot the electric scalar potential in close 2D space. The details in in this video and this one

The following in written in Matlab for quick prototype.

    % solve matrix
    clearvars
    % number of rows ans columns
    m = 21
    n = 41
    
    % m = 3
    % n = 3
    % initialise arrays
    U = zeros(m*n,m*n);
    
    dx = 0.01;
    dy = 0.01;
    
    for i = 1:m
        for j = 1:n
            index = (i-1)*n+j;
            U(index, index) = -2/dx^2 -2/dy^2;
            if j-1 > 0
                U(index,index-1) = 1/dx^2;
            end
    
            if j+1 <= n 
                U(index,index+1) = 1/dx^2;
            end
    
            if i-1 > 0
                U(index,index-3) = 1/dy^2;
            end
    
            if i+1 <= m
                U(index,index+3) = 1/dy^2;
            end
        end
    end
    
    % U
    
    V0 = 0.5
    % U*results =B
    B = zeros(m*n,1);
    B((m-1)*n+1:m*n, 1) = -V0/dy^2;
    
    
    % Solve Linear Equations in Matrix Form
    results = linsolve(U,B);
    
    final =zeros(m,n);
    
    for k = 1:m
        final(k,:) = results((k-1)*n+1:k*n);
    end    
    
    % plot final results
    f1 = figure('Name','Voltage surface plot','NumberTitle','off');
    [X,Y] = meshgrid(0:0.1:4, 0:0.1:2);
    s = surf(X,Y,final,'FaceAlpha',0.5)
    title("Potential V(x,y) for finite conducting box, " + "V_0=" + V0 + " volts")

    f2 = figure('Name','imagesc','NumberTitle','off');
        
    clims = [4 18];
    imagesc(X(1,:), Y(:,1), final)
    title("Potential V(x,y) for finite conducting box, " + "V_0=" + V0 + " volts")

Image from the author.

Image from me.