matlab code for a 2D steady state conduction problem using finite differencing method?

how can i get a matlab code for a 2D steady state conduction problem using finite differencing method? i want to check the heat flow at each of the nodes and edges and sum them to check if it is zero. i have got this code. how can i modify it to do what i want?

CODE: % Variable List: % T = Temperature (deg. Celsius) % T1 = Boundary condition temperature 1 (deg. Celsius) % T2 = Boundary condition temperature 2 (deg. Celsius) % Tinf = Ambient fluid temperature (deg. Celsius) % theta = Non-dimensionalized temperature difference = (T-T1)/(T2-T1) % Lx = Plate length in x-direction (m) % Ly = Plate length in y-direction (m) % AR = Aspect ratio of Ly / Lx to ensure dx = dy % h = Convection coefficient (W/m^2K) % k = Thermal conductivity (W/mK) % Bi = Finite-difference Biot number % Nx = Number of increments in x-direction % Ny = Number of increments in y-direction % dx = Increment size in x-direction (m) % dy = Increment size in y-direction (m) % dT = Temperature step between contours % tol = Maximum temperature difference for convergence (deg. Celsius) % pmax = Maximum number of iterations % Told = Stores temperature matrix for previous time step % diff = Maximum difference in T between iterations (deg. Celsius) % i = Current column % j = Current row % p = Current iteration % v = Sets temperature levels for contours % x = Create x-distance node locations % y = Create y-distance node locations % Nc = Number of contours for plot

clear, clc % Clear command window and workspace

Lx = .50; % Plate half-length in x-direction (m)

Ly = .20; % Plate length in y-direction (m)

Nx = 25; % Number of nodes in y-direction

Ny=25; % Number of nodes in x-direction

AR = Ly/Lx; % Aspect ratio of Ly /Lx to ensure dx = dy

%Ny = AR*Nx; % Number of increments in y-direction

dx = Lx/Nx; % Increment size in x-direction (m)

dy = Ly/Ny; % Increment size in y-direction (m)

T1 = 10; % BC temperature at end of fin (deg. Celsius)

T2 = 25; % BC temperature at base of fin (deg. Celsius)

Tinf = T2; % Ambient fluid temperature (deg. Celsius)

h = 30; % Convection coefficient (W/m^2K)

k = 10; % Thermal conductivity (W/m^2K)

Bi = h*dx/k; % Finite-difference Biot number

T = T1*ones(Nx+1,Ny+1); % Initialize T matrix to T1 everywhere

T(1:Nx+1,1) = T1; % Initialize base of fin to T1 BC

tol = 10^-2; % Max temp delta to converge (deg. Celsius) pmax = 10*10^2; % Maximum number of iterations

x = 0:dx:Lx; % Create x-distance node locations y = 0:dy:Ly; % Create y-distance node locations

for p = 1:pmax % Loop through iterations Told = T; % Store previous T array as Told for later

    for j = 2:Ny            % Loop through rows
        for i = 2:Nx        % Loop through interior columns

            % Calculates convection BC along left side
            if i == 2
                T(1,j) = (2*T(2,j)+T(1,j-1)+T(1,j+1)+2*Bi*Tinf)/(4+2*Bi); 
            end

            % Calculates interior node temperatures
            T(i,j) = .25*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1)); 

            % Calculates adiabatic BC along right side
            if i == Nx
                T(Nx+1,j) = .25*(2*T(Nx,j)+T(Nx+1,j-1)+T(Nx+1,j+1)); 
            end

        end
    end

    for i = 2:Nx  % Loop through interior columns at top row

        % Calculates top left corner, conv/conv BC's            
        if i == 2
            T(1,Ny+1) = (T(1,Ny)+T(2,Ny+1)+2*Bi*Tinf)/(2+2*Bi); 
        end

        % Calculates convection BC along top
        T(i,Ny+1) = (2*T(i,Ny)+T(i-1,Ny+1)+T(i+1,Ny+1)+2*Bi*Tinf)...
            /(4+2*Bi);

        % Calculates top right, conv/adiabatic BC's
        if i == Nx
            T(Nx+1,Ny+1) = (T(Nx+1,Ny)+T(Nx,Ny+1)+Bi*Tinf)/(2+Bi); 
        end

    end

    diff = max(max(abs(T - Told)));     % Max difference between iterations
    fprintf('Iter = %8.0f - Dif. = %10.6f deg. C\n', p, diff);
    if (diff < tol)     % Exit iteration loop because of convergence
        break
    end
end

fprintf('Number of iterations = \t %8.0f \n\n', p) % Print time steps

if (p == pmax) % Warn if number of iterations exceeds maximum disp('Warning: code did not converge') fprintf('\n') end

disp('Temperatures in block in deg. C = ') for j = Ny+1:-1:1 % Loop through each row in T fprintf('%7.1f', T(:,j)) % Print T for current row fprintf('\n') end

Nc = 50; % Number of contours for plot dT = (T2 - T1)/Nc; % Temperature step between contours v = T1:dT:T2; % Sets temperature levels for contours colormap(jet) % Sets colors used for contour plot contourf(x, y, T',v, 'LineStyle', 'none') colorbar % Adds a scale to the plot axis equal tight % Makes the axes have equal length title('Contour Plot of Temperature in deg. C') xlabel('x (m)') ylabel('y (m)') pause(0.001) % Pause between time steps to display graph figure surf(x,y, T') title('3D Projection of Temperature in deg. C')