How do i modify the code in order to make the boundary conditions variable (changing every iteration) for 1D transient heat transfer implicit finite difference method ?

clear all; 

% input data 
Kair=0.0242; %thermal conductivity of air 
Den_air=1.225; %density of air
Cp_air=1006.43; %specific heat of air
Kglass=0.96; % thermal conductivity of glass
Den_glass=2500; %density of glass
Cp_glass=840; %specific heat of glass
L=0.018; 
tf=20; %time limit 
Ta=23; %ambient temp (celcius)
Ts=40; %source heat (celcius)

%step 
Nt=101;
Nz=29;
t=linspace(0,tf,Nt);
z=linspace(0,L,Nz);
dt=t(2)-t(1);
dz=z(2)-z(1);

%Matrix prep
A=zeros(Nz,Nz);   %position
B=zeros(Nz,1);    %constant 
T=zeros(Nt,Nz);   %our goal

%fill T with IC 
%initial Condition 
T(1,:)=ones(1,Nz)*Ta; 

%BC-1
A(1,1)=1;
B(1)=Ts;

%BC-2
A(end,end)=1;
B(end)=Ta;

%for simplfy calculation
alfa_air=Kair/(Den_air*Cp_air)

alfa_air = 1.9629e-05

alfa_glass=Kglass/(Den_glass*Cp_glass)

alfa_glass = 4.5714e-07

var=dt./dz.^2; 

%tridiagonal matrix 
for j=2:Nt  %j=time, time start for 2
    for i=2:Nz-22 %i=position 
        s= alfa_glass*var; 
        A(i,i-1)=-s;
        A(i,i)=1+2*s; 
        A(i,i+1)=-s;
        B(i)=T(j-1,i); %j-1 means 2-1=1
    end 
    for i=8:Nz-6  %i=position 
        s= alfa_air*var; 
        A(i,i-1)=-s;
        A(i,i)=1+2*s; 
        A(i,i+1)=-s;
        B(i)=T(j-1,i); %j-1 means 2-1=1
    end    
    for i=24:Nz-1  %i=position 
        s= alfa_glass*var;
        A(i,i-1)=-s;
        A(i,i)=1+2*s; 
        A(i,i+1)=-s;
        B(i)=T(j-1,i); %j-1 means 2-1=1
    end 

T(j,:)=A\B;
figure (1)
plot(z,T(j,:), 'LineWidth',3)
title(['Temp at t=', num2str(t(j)),'sec' ])
axis ([0 L 0 Ts])
xlabel('Distance')
ylabel('Temperature')
pause (0.1)
end 

figure (2)
imagesc(z,t,T)
title ('Distribution Temperature')
xlabel('Distance')
ylabel('Time')
colormap jet
colorbar
grid on

the code above solves 1D transient heat transfer implicit form, but the boundary conditions are fixed. the attached pictures shows my boundary conditions where T1, T0, TNz, TNz-1 are variables they are chaning with time just like the interior nodes, the other values are all known like Tin, h, etc.. how do i modify the code based on these boundray conditions?

[Torsten]

Use "pdepe" with the internal nodes where the material properties change being part of the mesh.