# Problem with my Octave code (unsteady heat equation with FEM)

I want help with my Octave code regarding the unsteady heat equation.

My geometry and mesh are generated with FreeFEM++, so there is no problem with that (I tried it with the steady problem with no for loop and it worked just fine).

But when I introduce the for loop to iterate with the time step and try to visualize the results at multiple steps, the results remain the same for the entire loop with no change in figures or in Temperature data.

I tried varying the time step, the maximum number of iterations, the source term, and alpha (implicit/explicit Euler and Crank-Nicolson)

Here is the script :

% Visual aspects
set(0,'defaultAxesFontSize',20);
set(0,'defaultAxesLineWidth',2);

% Physical parameters
kappa=0.005;
amplitude_s=0.0; % 0.0
Thole=5;
Tborder=2;
Tcolormap=[2 5];

% Read the informations from the geometry

% Read matrices from the Finite Element Formulation
[n,m]=size(M);
%pcolor(M);
%
%
% Get the position of the Dirichlet masks
ii_Dir_hole=find(abs(BC_Dir_hole)>0);
ii_Dir_borders=find(abs(BC_Dir_borders)>0);
ii_int=find(BC_Dir==0);
%
% Build corresponding matrices
T=zeros(n,1);
alpha=0;
%alpha=1
%alpha=1/2
N=100;
delta_t=0.001;
A=1/delta_t*M+alpha*kappa*K;
Ared=A(ii_int,ii_int);
s=amplitude_s*ones(n,1);
s_new=s;
B=(1/delta_t*M-(1-alpha)*kappa*K)*T+M*(alpha*s_new+(1-alpha)*s);
% SOLVE
Td=zeros(n,1);
Td(ii_Dir_hole)=Thole;
Td(ii_Dir_borders)=Tborder;
%
for nt=1:N
C=B-A*Td;
Cred=C(ii_int);
Tred=Ared\Cred;
%
T=Td;
T(ii_int)=Tred;
%
if mod(nt,10)==0;
FIG=visu_2D(x,y,mesh_info,T,'none','off');
%FIG=visu_2D_visibility(x,y,mesh_info,T,'none','on');
caxis(Tcolormap);
%
savename=['PIC/PIC_' sprintf('%05d',not)];
print(FIG,'-r200','-djpg',savename);
close all;
end
endfor


Reading quickly through your code, it seems that you never update B. So the term $$T^n$$ of the previous step remains unchanged, hence you are effectively computing the same step over and over again...