1D diffusion equation of Heat Equation

I am a little confused about the output of the function t_crit, I have no ideas how to write this code, I think others are okay which is shown below.

and the question is shown like the graph below

function t_crit = function1(nx,dt,mFlag)
%% Solution to the 1D diffusion equation
% inputs:  - nx      Number of points in the domain
%          - dt      Temporal step size
%          - mFlag   the flag for band material
% outputs: - t_crit  how long the band stays above 43 degrees
% Define the parameters
r       = 0.03;             % Fixed Radius
L       = r*pi;             % Length of the domain
xstart  = -L;               % Start of computational domain (m)
xend    = xstart+2*L;       % End of computational domain (m)
T1      = 38;               % Initial temperature for most of the band
T2      = 50;               % Temperature in the square wave
phw      = pi/9;            % Pulse half-width, referring tp initial condition
nt = 101;
if mFlag == 0                            % Diffusivity for stainless/rubber m^2/s Tc - c^2Txx = 0 D = k/(sigma*rho)
    D = 14.4/(500*7900);        % diffusion of stanliess stee
elseif mFlag == 1
    D = 0.19/(2176*1230);  %diffusion of rubber
end

c       = sqrt(D) ;               % c coefficient in the heat equation

% Create the x grid
x       = linspace(xstart,xend,nx);
dx      = (xend-xstart)/(nx-1);

% Create the time grid 
timestart = 0;              % Starting time (0 seconds)
timeend = (nt-1)*dt+timestart;             % Temporal step size, delta t (s)
time      = linspace(timestart,timeend,nt);             % Time vector 
% alternatively: time = timestart :dt:timeend

% Fourier term parameters
nfs = 1000;                 % Number of terms in the Fourier Series

% Loop through time
for i = 1:nt                % For each time

    t = time(i);            % time in seconds

    % Initialise the solution vector
    T = zeros(1,nx);

    % Compute the solution at each point in domain
    for j = 1:nx
        % Compute steady-state solution at point
        A0 = 38 + 12/9;
        T(j) = A0;

        % Sum all terms in the homogeneous solution at point
        for n = 1:nfs
            An = 24/(n*pi)*sin(n*pi/9);         % Find the An at that n
            lambdan= (n*pi*c/L);     % Find the lambdan at that n
            T(j) = T(j) + An*exp(-lambdan^2*t);
        end
    end
end

% the input file of the function looks like 
% nx = 101;
% dt = 10;
% mFlag = 1; 
% t_crit = function1(nx,dt,mFlag)

[darova]

I think you should return matrix of size 101x1000 (every solution)