2-D Steady State Heat Conduction with Multiple Materials

Hello, I am trying to setup a Matlab code to solve a 2-D steady state heat conduction equation using the finite difference method. This is a picture of what I am trying to model:

This is the code I have written so far:

%initialize number of nodes used
Nx = 30;
Ny = 30;
Lx = 650;
Ly = 175;

%k and h values
k_a = 0.15;
k_b = 0.04;
k_c = 0.20;
k_d = 0.25;
h_a = 50;
h_b = 5;

%initial temperatures
T_upper = 323;
T_lower = 298;

%initialize function
T = zeros(Nx,Ny);

%increment values
dx = Lx/Nx;
dy = Ly/Ny;

%initialize boundary conditions
T(1,:) = T_upper; %sets bottom condition
T(Nx,:) = T_lower; %sets top condition

%calculate temperature distribution
flag = 1;
while flag == 1
    flag = 0;
    T_old = T;
    for i = 2:Nx-1
        for j = 2:Ny-1
            if i == 2
                %solves for the initial surface
                T(i,j) = ((-2*T_old(i,j+1) - T_old(i+1,j) - T_old(i-1,j)) -   ((2*h_a*dx*T_upper)/k_a))/(-2*(((h_a*dx)/k_a)+2));
                error = abs((T(i,j) - T_old(i,j))/T_old(i,j)); 
            elseif i>2 && i<12
                %solves  for block A
                T(i,j) = 0.25*(T_old(i-1,j) + T_old(i+1,j) + T_old(i,j-1) + T_old(i,j+1));
                error = abs((T(i,j) - T_old(i,j))/T_old(i,j));
            elseif i == 12
                %solves for surface of block B and C
                T(i,j) = ((-2*T_old(i,j+1) - T_old(i+1,j) - T_old(i-1,j)) -((2*h_a*dx*T(i,j-1))/k_b))/(-2*(((h_a*dx)/k_b)+2));
                error = abs((T(i,j) - T_old(i,j))/T_old(i,j));
            elseif i>12 && i<25
                %solves for block B and C
                T(i,j) = 0.25*(T_old(i-1,j) + T_old(i+1,j) + T_old(i,j-1) + T_old(i,j+1));
                error = abs((T(i,j) - T_old(i,j))/T_old(i,j));  
            elseif i>25 && i<Nx-1
                %solves for surface of block D
                T(i,j) = ((-2*T_old(i,j+1) - T_old(i+1,j) - T_old(i-1,j)) -((2*h_a*dx*T(i,j-1))/k_c))/(-2*(((h_a*dx)/k_c)+2));
                error = abs((T(i,j) - T_old(i,j))/T_old(i,j));  
            end
            if error > 0.01
                flag = 1;
            end
        end
    end
end

%plot temperature distribution
[x,y] = meshgrid(0:1/(Nx-1):1 , 0:1/(Ny-1):1);
contourf(x,y,T);
colorbar
colormap
title('Block A Temp Distribution')
xlabel('x')
ylabel('y')

What I am trying to do is model the above picture, but only using the boxed in version shown since it is repeating. The left and right sides will be adiabatic insulated. The top and bottom surfaces are exposed to the conditions shown.

I created this code initially for just a solid block and it runs correctly and gives me the correct plot. I am having trouble now updating that code with the new case of 4 different materials. I know some of my values are not in the correct units, but right now I am just trying to get the code to work and will fix the constants later.

Right now what I am trying to get the code to do is solve for first surface, then solve for block A according to the first surface, then solve for the interface of block A to block B/C, and so on and so forth.

I am also having trouble implementing how to differentiate block B and C since they have two different k values.

Can someone point me in the right direction?