2D Enclosure Problem Net Radiation Method

Hey guys. I have a problem. The code and the problem is given below. Does anyone help me about how to code the rest? I have the formulation for two-parallel surfaces but I cannot implement that into enclosure problem. If somebody help me to establish at least the pseudocode for establishing the linear system?

The view factor values are stored in "F" variable. Grid=2 means I divided each surface at 2.

clear;
clc;
x = [0 200 200 0 ];      % x-coordinates of end points
y = [0 0 100 100 ];      % y-coordinates of end points
x = [x 0]; y=[y 0];      
%%
% Calculating the number of enclosure surface
n_s = size(x,2)-1;
%%
% Calculating the angle and length
for i = 1:n_s
    %%
    % Lengths
    L(i,1) = sqrt((x(i+1)-x(i))^2+(y(i+1)-y(i))^2);
    %%
    % Angles
    thetac(i,1) = acosd((x(i+1)-x(i))/L(i));
    thetas(i,1) = asind((y(i+1)-y(i))/L(i));
end
%% Dicretization
% Nummber of grids for each surface
grid=2;
%%
% Evaluating the grid spacing for each surface
for i=1:n_s
    dL(i,1) = L(i)/grid;
end

%% Coordinates of Grid Points
for i=1:n_s
    xi(i,1)=x(i);
    yi(i,1)=y(i);
    
    for j=2:(grid+1)
        
        xi(i,j) = xi(i,j-1)+dL(i)*cosd(thetac(i));
        yi(i,j) = yi(i,j-1)+dL(i)*sind(thetas(i));
        
    end
end

%%
% Collecting the row values into a column vector
k=0;
for i=1:n_s
    for j=2:(grid+1)
        k=k+1;
        xx(k)=xi(i,j);
        yy(k)=yi(i,j);
    end
end
xx=[0 xx];
yy=[0 yy];
%% The view factor loop
for i=1:n_s*grid
    %%
    % The length of the originating surface
    dL_n(i) =sqrt(( xx(i+1)-xx(i))^2+(yy(i+1)-yy(i))^2); 
    %%
    for j = 1:n_s*grid
        if i==j
            F(i,j)=0;
        else
            %%
            % Calculating the diagonals
            diagonal_1= sqrt((xx(j+1)-xx(i+1))^2+(yy(j+1)-yy(i+1))^2);
            diagonal_2=sqrt((xx(j)-xx(i))^2+(yy(j)-yy(i))^2);
            %%
            % Calculating the Sides         
            side_1=sqrt((xx(j+1)-xx(i))^2+(yy(j+1)-yy(i))^2);
            side_2=sqrt((xx(j)-xx(i+1))^2+(yy(j)-yy(i+1))^2);
            %%
            % Evaluating the view factor    
            F(i,j)=((diagonal_1+diagonal_2)-(side_1+side_2))/(2*dL_n(i));
        end
    end
end
%% Properties and Iteration
% Properties
T_1 = 1000;
T_3 = 500;
E_1 = 0.7;
E_2 = 1;
E_3 = 0.3;
E_4 = 1;
Q_1 = 0;
Q_4 = 0;
sigma = 5.670374419e-8;
% Matrices
A = zeros(n_s*grid,n_s*grid);
x = zeros(n_s*grid,1);
b = zeros(n_s*grid,1);

%% Iteration