Laplace equation using central difference

Hi im trying to program a code that solves laplace equation using central difference. I already set up the boundary condition that corresponds to b matrix in the code, however im having hard time coding the A matrix that contains all the coefficients that looks like following. The A and b are found in the for loop in the code.

close all
h = 1.25; %Try different values & find the order of accuracy
% Geometry & boundary conditions
xmin = 0.0; xmax = 40.0;
ymin = 0.0; ymax = 40.0;
% Boundary conditions (In this problem, all boundaries are isothermal)
Tl = 75; %left
Tr = 50; %right
Tb = 0; %bottom
Tt = 100; %top
% Mesh (h is given as input)
nx = (xmax-xmin)/h - 1; %number of interior nodes in x-direction
ny = (ymax-ymin)/h - 1; %number of interior nodes in y-direction
N = nx*ny; %total number of interior nodes
% Allocate memory for A and b (sparse matrix & vector).
A = sparse(N,N); b = sparse(N,1);
% A=zeros(N,N);
% b=zeros(N,1);
% Loop through all the interior nodes, and fill A & b.
for k=1:N
    A(k,k) = 4; % the diagonal element, corresponding to Tij, is always 4
    %In the following, we look at the four neighbours of (i,j) that are
    %involved in the 5-point formula
    i = mod(k-1,nx)+1; j = (k-i)/nx+1; %get (i,j) from k.
disp(i)
disp(j)
    if i==1
        b(k,1)=Tl;
    elseif i==nx
        b(k,1)=Tr;
        
    elseif j==1
        b(k,1)=Tb;
    elseif j==ny
        b(k,1)=Tt;
    end
    


end

% Solve Ax = b for nodal temperature values
T = full(A\b);

% Output Temperature at a sensor location (10 cm, 10 cm)
i = (0.25*(xmax-xmin))/h;
j = (0.25*(ymax-ymin))/h;
Tsensor = T(nx*(j-1)+i);
fprintf('h = %e, T(%d,%d) = %e.\n', h, i, j, Tsensor);

% Visualization
figure
[X, Y] = meshgrid(xmin:h:xmax, ymin:h:ymax);
Tvis = zeros(size(X)); %just for visualization
Tvis(:,1) = Tl; Tvis(:,nx+2) = Tr; Tvis(1,:) = Tb; Tvis(ny+2,:) = Tt;
for j = 2:nx+1
  for i = 2:ny+1
      Tvis(i,j) = T((i-2)*nx + j-1);
  end
end 
h = surf(X,Y,Tvis);
colormap('hot');
view(0,90);
shading interp
%set(h,'edgecolor','k');
cb = colorbar; cb.Label.String = 'Temperature (C)';
xlabel('X (cm)')
ylabel('Y (cm)')
zlabel('Temperature (C)')
daspect([1 1 1]);

[darova]

Lokk into this script. It should be helpfull

[Faraz Vossoughian]

Hi darobva, my question is set up differently, im really having a hard time codeing the coefficient matrix (A). could you please help

for k=1:N
    A(k,k) = 4; % the diagonal element, corresponding to Tij, is always 4
%     In the following, we look at the four neighbours of (i,j) that are
%     involved in the 5-point formula
    i = mod(k-1,nx)+1; j = (k-i)/nx+1; %get (i,j) from k.
    disp('ass')


    if i==1

        b(k,1)=Tl;
   
    elseif i==nx
        b(k,1)=Tr;
        
    elseif j==1
        b(k,1)=Tb;
    elseif j==ny
        b(k,1)=Tt;
    end
% if j<ny-2
%    r=3+(j-1)
%     A(i,nx-r)=-1
% end


%
if i>1 && j>1
if i==j

            A(j-1,i)=-1;
            A(j,i+1)=-1;
            A(j,i-1)=-1;
            A(j+1,i)=-1;
  
end
end
    
end

[darova]

Please look into the script i attached. Especially at these lines

for i = 2:m-1
    for j = 2:n-1
        ii = i + (j-1)*m;
        A(ii,ii-1) = ci;    % U(i-1,j) 
        A(ii,ii+1) = ci;    % U(i+1,j) 
        A(ii,ii)   = cij;   % U(i,j) 
        A(ii,ii-m) = cj;    % U(i,j-1) 
        A(ii,ii+m) = cj;    % U(i,j+1) 
    end
end