adding matrices into an element of another matrices

I have a for loop running representing time. I calcutale a value Tnew in the for loop the value then hase to be set back to Told when running the for loop again but then the values in Told becomes disregarded/ forgotten

thus I want to store the value into a matrics where each element represents a time step, I do not want to mesh the matrices a want to keep them seperate.

L = 4;
rg = 0.2;
total_time = 24*3600;
g = 9.81 ;

nx = 101;
nr = 51;
nt = 20;

dx = L/(nx-1);
dr = rg/(nr-1);
dt = total_time/(nt-1)

Told = ones(nx,nr)*(15 + 273.15);  % temperature at t = 0
Told(:,nr) = (65 + 273.15)
Tnew = zeros(nx,nr);
Ttime = cat(nt)
% m = 1;
% n = 1;
% t = 0;

r = 0:dr:rg
x = 0:dx:L
ts = 1

for t = 0:dt:total_time
    
    Tcomp = 65 + 273.15;
    Tamb = (12*sin((pi*t/12)-(2*pi/3)) + 8) + 273.15;
    pair = 101325/(287.08*Tamb);
    cpair = 1.045356e3 - (3.161783e-1)*Tamb + (7.083814e-4)*(Tamb)^2 - (2.705209e-7)*(Tamb)^3;
    uair = 2.287973e-6 + (6.259793e-8)*Tamb - (3.131956e-11)*(Tamb)^2 + (8.15038e-15)*(Tamb)^3;
    kair = -4.937787e-4 +(1.018087e-4)*Tamb - (4.627937e-8)*(Tamb)^2 +(1.250603e-11)*(Tamb)^3;
    vair = uair/pair;
    Pr = (uair*cpair)/kair;
    
    for n = 1:nr
        
        
        Ao = 2*pi*r(n)*dx + pi*dr*dx;
        Ai = 2*pi*r(n)*dx - pi*dr*dx;
        As = 2*pi*r(n)*dr;
        V_element = 2*pi*r(n)*dr*dx;
        
        
        
        for m = 1:nx
            
            kw = -0.614255 + (6.9962e-3)*Told(m,n) - (1.01075e-5)*(Told(m,n))^2 + (4.74737e-12)*(Told(m,n))^4;
            pw = 1/(1.49343e-3 - (3.7164e-6)*(Told(m,n)) + (7.09782e-9)*(Told(m,n))^2 - (1.90321e-20)*(Told(m,n))^6);
            cpw = 8.15599e3 - (2.80627e1)*Told(m,n) + (5.11283e-2)*(Told(m,n))^2 - (2.175882e-13)*(Told(m,n))^6;
            
            %Calculations for h -> convection coeficcient
            
            B = -(1/pw)*((pair - pw)/(Tamb - Told(m,n)));
            Lc = 2*rg;
            
            Ra = (g*B*(Tamb-Told(m,n))*(Lc^3)*Pr)/(vair^2);
            
            Nu = (0.6 + (0.387*(Ra^(1/6)))/(1 + (0.559/Pr)^(9/16))^(8/27))^2;
            
            h = (Nu*kair)/Lc;
            
            
            
         
            Tcomp = 65 + 273.15;
            
          
            
            Tnew(m,n) = ((Q_LHS + Q_RHS + Q_I + Q_O + Q_loss)*dt)/(div_by) + Told(m,n);
            
            
            
        end
        
        m = 1;
        
    end
    %     t = t + dt
    n = 1;
    Told = Tnew;                             %THIS IS WHERE I WOULD LIKE TO "SAVE" THE VALUE
 Ttime(ts) = Tnew;
   ts = ts+1; 
   
end