I get an incorrect T value. Could you help me with this please?

function [T,aliveFlag] = tutorial3(T1,timeend,nt)
% Description: This file implements the Fourier Series solution to the heat equation.
%
% Input parameters:
% - T1        - temperature at r=0
% - timeend   - time of simulation in seconds
% - nt        - number of points in time vector
%
% Output parameters:
% - T         - temperature profile along persons head
% - aliveFlag - 1 if person is alive, 0 if the person is dead
% Spatial (r) parameters
L      = 0.12;                      % Size of the domain
rstart = 0;                         % Start of computational domain (m)
rend   = L;                         % End of computational domain (m)
nr     = 101;                       % Number of grid points along the radius vector
dr     = (rend-rstart)/(nr-1);      % Spatial step size, delta r (m)
r      = linspace(rstart,rend,nr);            % Vector of grid points in r
% Temporal (t) parameters 
timestart = 0;% Starting time (0 seconds)
timeend = 60;
dt        = (timeend-timestart)/(nt-1);  % Temporal step size, delta t (s)
time      = linspace(timestart,timeend,nt);        % Vector of times
% Phyiscal parameters
T1    = -41;                        % Temperature at r=0
T2    = 37;                         % Temperature at r=L
A     = T2-T1;                      % Amplitude of the saw tooth wave
k     = 0.527;                      % Thermal conductivity (W/m-K)
rho   = 1000;                       % Density (kg/m^3)
gamma = 3600;  % Heat capacity (J/kg-
c     = sqrt(k/(rho*gamma));     % Heat equation constant   c^2=k/(rho*sigma)              
% Calculate B coefficients using a for loop
nfs = 1000;               % Number of fourier terms               
B = zeros(1,nfs);         % Initialise B vector   
lambda=zeros(1,nfs);      % Initialise lambda vector
for n = 1:nfs;           
    B(n)=156/(n*pi);                   % Calculate B coefficients
    lambda(n)=(n*pi*c)/L;
end
%% Solve for er series
% Loop thT using three for loops in time, space, and Fourirough time
for i = 1:nt               % For each time
    t = time(i);           % time in seconds
    % Vector of zeros to initialise the Fourier series solution.
    % This should be re-initialised at each new time step.
    T = zeros(1,nr);
    % Loop through space
    for j = 1:nr;          % For each grid point        
        T(j) = (78*r(j)/L)-41;         % Add the steady state solution 
        % Loop through the Fourier series
        for n = 1:nfs;
            % Calculate series sum for T at r(j) and t
         T(j)    = B(n)*sin(n*pi*r(j)/L)*exp(-lambda(n)^2*t)+T(j);             
        end
    end
end

    
aliveFlag=mean(T)>22;
disp(T);
disp(aliveFlag);

end

[KALYAN ACHARJYA]

Pass the inputs arguments tutorial3(T1,timeend,nt) and get the T, as defined output argument.

[Siratul Srotoshwini]

I did pass the input arguments but I still get an incorrect value. Could you please check if there is sth wrong with my code or not?