Improving Code (bisection method)

Hello everyone. How can I modify this script so that TL is no longer constant and v_avg no longer takes only 3 values, but instead has variable inputs over time taken from an Excel sheet?

%% Working code
clear;
clc;
global  M R_gas k D cp_G mu_G cpL h_L rho_liq DeltaH_evap H_c TL A_ant B_ant C_ant NumbersSpatialDiscretisation t_riempimento t_riposo t_svuotamento cycle_duration
%% inputs 
T_t = 328;  %268,283,288,298,328
%% Parametri del serbatoio e costanti fisiche
H = 15; % Altezza del serbatoio (m)
diametro = H; % Diametro del serbatoio (m)
V_serbatoio = pi * (diametro / 2)^2 * (H-1); % Volume del serbatoio (m^3)
P = 1e5; % Pressione in Pa
%% Costanti Diesel
M = 0.14; % Massa molare del Diesel in kg/mol
R_gas = 8.314; % Costante dei gas perfetti in J/(mol*K)
k = 0.025; % Conducibilità termica dell'aria in W/m·K
D = 10^(-5); % Diffusività in m²/s
cp_G = 2100; % Calore specifico in J/kg·K
mu_G = 1.85e-5; % Viscosità dinamica dell'aria (Pa·s)
cpL = 2200; % Calore specifico Diesel liquido in J/kg·K
h_L = 12; % Coefficiente scambio termico del liquido W/m²·K
rho_liq = 835; % Densità del Diesel liquido a 298 K in kg/m³
DeltaH_evap = 375000*M; % Calore latente di evaporazione in J/mol
H_c = 0; % Costante di Henry NON UTILIZZATA
TL = 283; % Temperatura liquido in K
%% Costanti Antoine (per tensione di vapore)
A_ant = 12.101;
B_ant = 8907;
C_ant =0;
%% Parametri discretizzazione
NumbersSpatialDiscretisation = 70; % Numero di nodi spaziali
%% Condizioni iniziali
T_init = T_t * ones(NumbersSpatialDiscretisation, 1); % Temperatura iniziale uniforme (K)
C_init = 0.16*ones(NumbersSpatialDiscretisation,1);
%C_init = linspace(0.02, 0.01, NumbersSpatialDiscretisation)'; % Gradiente iniziale della concentrazione
y0 = [C_init; T_init]; % Condizioni iniziali combinate
%% Parametri per i cicli di riempimento, riposo e svuotamento
t_riempimento = 24 * 3600; % 24 ore
t_riposo = 1 * 24 * 3600; % 5 giorni di riposo
t_svuotamento = 24 * 3600; % 24 ore di svuotamento
cycle_duration = t_riempimento + t_riposo + t_svuotamento;
transition_time = 3600; % Tempo di transizione (1 ora)
%% Funzione per calcolare v_avg
v_avg_func = @(t) v_cicli(t, V_serbatoio, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time);
%% Definizione dei parametri extra da passare
params = {NumbersSpatialDiscretisation, D, k, cp_G, mu_G, diametro, v_avg_func, P, M, R_gas, h_L, rho_liq, cpL, DeltaH_evap, A_ant, B_ant, C_ant, H, H_c,V_serbatoio, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time,TL,T_t};
%% Matrice di massa (per condizioni al contorno)
mass_matrix = eye(2*NumbersSpatialDiscretisation);
% initial point for C in space mesh
mass_matrix(1,1) = 0;
% final point for C in space mesh
mass_matrix(NumbersSpatialDiscretisation,NumbersSpatialDiscretisation) = 0;
% initial point for T in space mesh
mass_matrix(NumbersSpatialDiscretisation+1,NumbersSpatialDiscretisation+1) = 0;
% final point for T in space mesh
mass_matrix(2*NumbersSpatialDiscretisation,2*NumbersSpatialDiscretisation) = 0;
t_span = [0,cycle_duration*30]
%% Risoluzione del sistema tramite ode15s per ogni intervallo temporale
Output = OdeCompiler(t_span, params{:},mass_matrix,y0);
t = Output(:,1);
u_sol = Output(2:end,:);
%% FCN with time related jumps
function out = OdeCompiler(t_span, ~, ~, ~, ~, ~, diametro, ~, P, ~, ~, ~, ~, ~, ~, ~, ~, ~, H, ~,V_serbatoio, ~, ~, ~, ~, ~,~,T_t,~,~)
global  M R_gas k D cp_G mu_G cpL h_L rho_liq DeltaH_evap H_c TL A_ant B_ant C_ant NumbersSpatialDiscretisation t_riempimento t_riposo t_svuotamento cycle_duration
%% Parametri discretizzazione
y = linspace(0,1,NumbersSpatialDiscretisation);
%% Condizioni iniziali
T_init = TL * ones(NumbersSpatialDiscretisation, 1); % Temperatura iniziale uniforme (K)
T_init(NumbersSpatialDiscretisation) = T_t;
C_init = 0.01*ones(NumbersSpatialDiscretisation,1);
%C_init = linspace(0.02, 0.01, NumbersSpatialDiscretisation)'; % Gradiente iniziale della concentrazione
y0 = [C_init; T_init]; % Condizioni iniziali combinate
%% Parametri per i cicli di riempimento, riposo e svuotamento
t_riempimento = 12 * 3600; % 12 ore
t_riposo = 10 * 24 * 3600; % 7 giorni di riposo
t_svuotamento = 20*24 * 3600; % 24 ore di svuotamento
cycle_duration = t_riempimento + t_riposo + t_svuotamento+t_riposo/30;
transition_time = 3600; % Tempo di transizione (1 ora)
%% Funzione per calcolare v_avg
v_avg_func = @(t) v_cicli(t, V_serbatoio, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time);
%% Definizione dei parametri extra da passare
params = {NumbersSpatialDiscretisation, D, k, cp_G, mu_G, diametro, v_avg_func, P, M, R_gas, h_L, rho_liq, cpL, DeltaH_evap, A_ant, B_ant, C_ant, H, H_c,V_serbatoio, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time,TL,T_t};
%% Matrice di massa (per condizioni al contorno)
mass_matrix = eye(2*NumbersSpatialDiscretisation);
% initial point for C in space mesh
mass_matrix(1,1) = 0;
% final point for C in space mesh
mass_matrix(NumbersSpatialDiscretisation,NumbersSpatialDiscretisation) = 0;
% initial point for T in space mesh
mass_matrix(NumbersSpatialDiscretisation+1,NumbersSpatialDiscretisation+1) = 0;
% final point for T in space mesh
mass_matrix(2*NumbersSpatialDiscretisation,2*NumbersSpatialDiscretisation) = 0;
%compilazione del vettore tempo dividendolo in range con velocità del
%fluido costante per evitare cambi improvvisi durente la risoluzione ODE
NumberOfCycles = ceil(t_span(2)/cycle_duration);
Time_vector = zeros(NumberOfCycles*3+1,1);
for i = 0: (NumberOfCycles-1)
    Time_vector(i*3+2) = cycle_duration*i + t_riempimento;
    Time_vector(i*3+3) = cycle_duration*i + t_riempimento + t_riposo;
    Time_vector(i*3+4) = cycle_duration*i + t_riempimento + t_riposo +  t_svuotamento;
end
v_riempimento = V_serbatoio / (pi * (diametro / 2)^2 * t_riempimento);
v_svuotamento = -0.3*V_serbatoio / (pi * (diametro / 2)^2 * t_svuotamento);
timeSolution_vector = zeros(1,1);
Solution_vector = zeros(2,NumbersSpatialDiscretisation);
% calcolo velocità del fluido
for i = 2 : length(Time_vector)
    cycle_time = mod(Time_vector(i), cycle_duration+0.001);
    if cycle_time <= t_riempimento
        v_avg = v_riempimento;
    elseif cycle_time <= (t_riempimento + t_riposo)
        v_avg = 0;
    else
        v_avg = v_svuotamento;
    end
    % parametri ODE
    params = {NumbersSpatialDiscretisation, D, k, cp_G, mu_G, diametro, v_avg_func, P, M, R_gas, h_L, rho_liq, cpL, DeltaH_evap, A_ant, B_ant, C_ant, H, H_c,V_serbatoio, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time,TL,T_t,v_avg};
    tic 
    % opzioni ODE
    opts = odeset('RelTol',1e-10, 'AbsTol',1e-10, 'Mass', mass_matrix,'MassSingular','yes'); % Opzioni del solver
    if i > 2
        disp(i)
        pde_fun(t(end), u_sol(end,:)', params{:});
    else
    end
    % solutioni ODE
    [t,u_sol] = ode15s(@(t, u) pde_fun(t, u, params{:}), [Time_vector(i-1),Time_vector(i)], y0, opts);
    toc
    if i == 2
        Solution_vector = u_sol;
        timeSolution_vector = t;
    else
        Solution_vector = [Solution_vector; u_sol];
        timeSolution_vector = [timeSolution_vector;t];
    end
    y0 = u_sol(end,:)';
    
    pde_fun(t(1), u_sol(1,:)', params{:});
end
out = [timeSolution_vector,Solution_vector];
figure(1)
plot(timeSolution_vector/3600,Solution_vector(:,NumbersSpatialDiscretisation))
title('Andamento Ctop nel tempo')
xlabel('tempo [h]')
ylabel('Ctop [mol/m3]')
figure(2)
plot(y,Solution_vector(2000,1:NumbersSpatialDiscretisation))
title('Andamento C nello spazio')
xlabel('y[-]')
ylabel('C [mol/m3]')
figure(3)
plot(y,Solution_vector(2000,NumbersSpatialDiscretisation+1:2*NumbersSpatialDiscretisation))
title('Andamento T nello spazio')
xlabel('y[-]')
ylabel('T [K]')
figure(4)
surf(y,timeSolution_vector/3600,Solution_vector(:,1:NumbersSpatialDiscretisation))
title('grafico 3D C')
xlabel('y[-]')
ylabel('tempo [h]')
zlabel('C [mol/m3]')
figure(5)
surf(y,timeSolution_vector/3600,Solution_vector(:,NumbersSpatialDiscretisation+1:2*NumbersSpatialDiscretisation))
title('grafico 3D T')
xlabel('y[-]')
ylabel('tempo [h]')
zlabel('T [K]')
pde_fun(timeSolution_vector(1), Solution_vector(1,:),params{:});
end
%% Funzione PDE con discrezione spaziale
function Output_array = pde_fun(t, u, NumbersSpatialDiscretisation, D, k_G, cp_G, mu_G, diametro, ~, P, M, R_gas, ~, ~, ~, ~, A_ant, B_ant, ~, H, ~,~, t_riempimento, t_riposo, ~, cycle_duration, ~,TL,T_t,v_avg)
%adattamento nel tempo per fvorire stabilità solver durante riempimento e
%svuotamento
number_cycles = floor(t/cycle_duration);
if v_avg < 0
    t = t - (t_riempimento+t_riposo)*number_cycles;
elseif v_avg > 0
    t = t - cycle_duration*(number_cycles);
else
end
% Ricezione input
Concentration_array = u(1:NumbersSpatialDiscretisation);
Temperature_array = u((NumbersSpatialDiscretisation+1):(2*NumbersSpatialDiscretisation));
% Discretizzazione derivate spaziali
dTdy = zeros(NumbersSpatialDiscretisation,1);
d2Tdy2 = zeros(NumbersSpatialDiscretisation,1);
dCdy = zeros(NumbersSpatialDiscretisation,1);
d2Cdy2 = zeros(NumbersSpatialDiscretisation,1);
dy = 1/NumbersSpatialDiscretisation;
y = linspace(0,1,NumbersSpatialDiscretisation);
for i = 2:(NumbersSpatialDiscretisation-1)
    % Calcolo della derivata della concentrazione con il flux limiter
    rC = (Concentration_array(i) - Concentration_array(i-1)) / (Concentration_array(i+1) - Concentration_array(i));
    limiterC = flux_limiter(rC, 'minmod');  % Usa il limiter scelto
    
    dCdy(i) = limiterC * (Concentration_array(i+1) - Concentration_array(i)) / dy;  % Differenza limitata
    
    % Calcolo della derivata della temperatura con il flux limiter
    rT = (Temperature_array(i) - Temperature_array(i-1)) / (Temperature_array(i+1) - Temperature_array(i));
    limiterT = flux_limiter(rT, 'minmod');  % Usa il limiter scelto
    
    dTdy(i) = limiterT * (Temperature_array(i+1) - Temperature_array(i)) / dy;  % Differenza limitata
    
    % derivate seconde solo con central
    d2Tdy2(i) = (Temperature_array(i+1) - 2 * Temperature_array(i) + Temperature_array(i-1)) / dy^2;
    d2Cdy2(i) = (Concentration_array(i+1) - 2 * Concentration_array(i) + Concentration_array(i-1)) / dy^2;
end
dCdy(1) = (Concentration_array(2) - Concentration_array(1)) / (dy);
dTdy(NumbersSpatialDiscretisation) = (Temperature_array(NumbersSpatialDiscretisation) - Temperature_array(NumbersSpatialDiscretisation-1)) / (dy);
dCdy(NumbersSpatialDiscretisation) = (Concentration_array(NumbersSpatialDiscretisation) - Concentration_array(NumbersSpatialDiscretisation-1)) / (dy);
% Calcolo C saturazione con Raoult assumendo gas perfetti
SaturationConcentration = 6894.76*(exp(A_ant - B_ant / (1.8*Temperature_array(1))))/P* P/(R_gas*Temperature_array(1));
Interface_gasDensity=P/(R_gas*Temperature_array(1));
% Definisco v*
v_star = v_avg-(D/Interface_gasDensity).*dCdy(1);
% inizializzazione di un calcolo itarativo per la quantificazione dei
% coefficienti di trasporto
Speed_array = zeros(NumbersSpatialDiscretisation,1);
Gas_density = Speed_array;
Re = Speed_array;
Sc = Speed_array;
Pe_m = Speed_array;
Pe_t = Speed_array;
E_m = Speed_array;
E_t = Speed_array;
alfa = Speed_array;
Gas_density(1) = Interface_gasDensity;
%metodo della bisezione per trovare v
SpeedPreviousIteration = v_star;
LowerLimit = SpeedPreviousIteration/5;
Re(1) =LowerLimit .* diametro.* Gas_density(1) ./ mu_G; % Numero di Reynolds
Re(1) = abs(Re(1));
Pr = mu_G * cp_G / k_G; % Numero di Prandtl
Sc(1) = mu_G / (Gas_density(1)*M * D); % Numero di Schmidt
Pe_m(1) = (1 / (Re(1) .* Sc(1)) + (Re(1) .* Sc(1)) / 192); % Numero di Peclet materiale
Pe_t(1)= (1 / (Re(1) * Pr) + (Re(1) * Pr) / 192); % Numero di Peclet termico
E_m(1)=Pe_m(1)*D; %Coeff. dispersione materiale m^2/s
E_t(1)=Pe_t(1)*k_G./(Gas_density(1)*M*cp_G); %Coeff. dispersione termica m^2/s
alfa_LL = E_t(1)/Interface_gasDensity/cp_G;
Delta_Speed_LL = LowerLimit - (v_star+alfa_LL./Temperature_array(1).*(dTdy(i)-dTdy(end)));
UpperLimit = SpeedPreviousIteration*5;
Re(1) =UpperLimit .* diametro.* Gas_density(1) ./ mu_G; % Numero di Reynolds
Re(1) = abs(Re(1));
Pr = mu_G * cp_G / k_G; % Numero di Prandtl
Sc(1) = mu_G / (Gas_density(1)*M * D); % Numero di Schmidt
Pe_m(1) = (1 / (Re(1) .* Sc(1)) + (Re(1) .* Sc(1)) / 192); % Numero di Peclet materiale
Pe_t(1)= (1 / (Re(1) * Pr) + (Re(1) * Pr) / 192); % Numero di Peclet termico
E_m(1)=Pe_m(1)*D; %Coeff. dispersione materiale m^2/s
E_t(1)=Pe_t(1)*k_G./(Gas_density(1)*M*cp_G); %Coeff. dispersione termica m^2/s
itCount = 0;
Speed_array(1) = LowerLimit + (UpperLimit - LowerLimit)/2;
Re(1) =Speed_array(1) .* diametro.* Gas_density(1) ./ mu_G; % Numero di Reynolds
Re(1) = abs(Re(1));
Pr = mu_G * cp_G / k_G; % Numero di Prandtl
Sc(1) = mu_G / (Gas_density(1)*M * D); % Numero di Schmidt
Pe_m(1) = (1 / (Re(1) .* Sc(1)) + (Re(1) .* Sc(1)) / 192); % Numero di Peclet materiale
Pe_t(1)= (1 / (Re(1) * Pr) + (Re(1) * Pr) / 192); % Numero di Peclet termico
E_m(1)=Pe_m(1)*D; %Coeff. dispersione materiale m^2/s
E_t(1)=Pe_t(1)*k_G./(Gas_density(1)*M*cp_G); %Coeff. dispersione termica m^2/s
alfa(1) = E_t(1)/Interface_gasDensity/cp_G;
delta_Speed = Speed_array(1) - (v_star+alfa(i)./Temperature_array(1).*(dTdy(i)-dTdy(end)));
% solving the equations to find the true vapour speed at the 
% vapour-liquid interface using the bisection method
while abs(delta_Speed)/abs(Speed_array(1)) > 1e-3 && itCount <100
    itCount = itCount + 1;
    Speed_array(1) = LowerLimit + (UpperLimit - LowerLimit)/2;
    % SpeedPreviousIteration = Speed_array(1);
    Re(1) =Speed_array(1) .* diametro.* Gas_density(1) ./ mu_G; % Numero di Reynolds
    Re(1) = abs(Re(1));
    Pr = mu_G * cp_G / k_G; % Numero di Prandtl
    Sc(1) = mu_G / (Gas_density(1)*M * D); % Numero di Schmidt
    Pe_m(1) = (1 / (Re(1) .* Sc(1)) + (Re(1) .* Sc(1)) / 192); % Numero di Peclet materiale
    Pe_t(1)= (1 / (Re(1) * Pr) + (Re(1) * Pr) / 192); % Numero di Peclet termico
    E_m(1)=Pe_m(1)*D; %Coeff. dispersione materiale m^2/s
    E_t(1)=Pe_t(1)*k_G./(Gas_density(1)*M*cp_G); %Coeff. dispersione termica m^2/s
    alfa(1) = E_t(1)/Interface_gasDensity/cp_G;
    delta_Speed = Speed_array(1) - (v_star+alfa(i)./Temperature_array(1).*(dTdy(i)-dTdy(end)));
    if delta_Speed*Delta_Speed_LL < 0
        UpperLimit = Speed_array(1);
        
    else
        LowerLimit = Speed_array(1);
        Delta_Speed_LL = delta_Speed;
    end
end
for i = 2:(NumbersSpatialDiscretisation)
    % initialisation of the bisection method for finding the true speed at
    % a given point in space
    % UL calcs
    Speed_UL = Speed_array(1)*5;
    Gas_density(i) = P/(R_gas*(TL + T_t)/2);
    Re(i) =Speed_UL .* diametro.* Gas_density(i) ./ mu_G; % Numero di Reynolds
    Re(i) = abs(Re(i));
    Pr = mu_G * cp_G / k_G; % Numero di Prandtl
    Sc(i) = mu_G / (Gas_density(i)*M * D); % Numero di Schmidt
    Pe_m(i) = (1 / (Re(i) .* Sc(i)) + (Re(i) .* Sc(i)) / 192); % Numero di Peclet materiale
    Pe_t(i)= (1 / (Re(i) * Pr) + (Re(i) * Pr) / 192); % Numero di Peclet termico
    E_m(i)=Pe_m(i)*D; %Coeff. dispersione materiale m^2/s
    E_t(i)=Pe_t(i)*k_G./(Gas_density(1)*M*cp_G); %Coeff. dispersione termica m^2/s
    alfa(i) = E_t(i)/Interface_gasDensity/cp_G;
    
    % LL calcs
    Speed_LL = Speed_array(1)/5;
    Gas_density(i) = P/(R_gas*Temperature_array(i));
    Re(i) =Speed_LL .* diametro.* Gas_density(i) ./ mu_G; % Numero di Reynolds
    Re(i) = abs(Re(i));
    Pr = mu_G * cp_G / k_G; % Numero di Prandtl
    Sc(i) = mu_G / (Gas_density(i)*M * D); % Numero di Schmidt
    Pe_m(i) = (1 / (Re(i) .* Sc(i)) + (Re(i) .* Sc(i)) / 192); % Numero di Peclet materiale
    Pe_t(i)= (1 / (Re(i) * Pr) + (Re(i) * Pr) / 192); % Numero di Peclet termico
    E_m(i)=Pe_m(i)*D; %Coeff. dispersione materiale m^2/s
    E_t(i)=Pe_t(i)*k_G./(Gas_density(1)*M*cp_G); %Coeff. dispersione termica m^2/s
    alfa(i) = E_t(i)/Interface_gasDensity/cp_G;
    delta_Speed_LL = Speed_LL - (v_star+alfa(i)./Temperature_array(i).*(dTdy(i)-dTdy(end)));
    %% calcs for i 
    Speed_array(i) = Speed_LL + (Speed_UL - Speed_LL)/2;
    Gas_density(i) = P/(R_gas*Temperature_array(i));
    Re(i) =Speed_array(i) .* diametro.* Gas_density(i) ./ mu_G; % Numero di Reynolds
    Re(i) = abs(Re(i));
    Pr = mu_G * cp_G / k_G; % Numero di Prandtl
    Sc(i) = mu_G / (Gas_density(i)*M * D); % Numero di Schmidt
    Pe_m(i) = (1 / (Re(i) .* Sc(i)) + (Re(i) .* Sc(i)) / 192); % Numero di Peclet materiale
    Pe_t(i)= (1 / (Re(i) * Pr) + (Re(i) * Pr) / 192); % Numero di Peclet termico
    E_m(i)=Pe_m(i)*D; %Coeff. dispersione materiale m^2/s
    E_t(i)=Pe_t(i)*k_G./(Gas_density(1)*M*cp_G); %Coeff. dispersione termica m^2/s
    alfa(i) = E_t(i)/Interface_gasDensity/cp_G;
    delta_Speed = Speed_array(i) - (v_star+alfa(i)./Temperature_array(i).*(dTdy(i)-dTdy(end)));
    itCount = 0;
    % solving the equations to find the true vapour speed at 
    % a given point in the mesh using the bisection method
    while abs(delta_Speed/Speed_array(1)) > 1e-3 && itCount <100
        itCount = itCount + 1;
        Speed_array(i) = Speed_LL + (Speed_UL - Speed_LL)/2;
        Gas_density(i) = P/(R_gas*Temperature_array(i));
        Re(i) =Speed_array(i) .* diametro.* Gas_density(i) ./ mu_G; % Numero di Reynolds
        Re(i) = abs(Re(i));
        Pr = mu_G * cp_G / k_G; % Numero di Prandtl
        Sc(i) = mu_G / (Gas_density(i)*M * D); % Numero di Schmidt
        Pe_m(i) = (1 / (Re(i) .* Sc(i)) + (Re(i) .* Sc(i)) / 192); % Numero di Peclet materiale
        Pe_t(i)= (1 / (Re(i) * Pr) + (Re(i) * Pr) / 192); % Numero di Peclet termico
        E_m(i)=Pe_m(i)*D; %Coeff. dispersione materiale m^2/s
        E_t(i)=Pe_t(i)*k_G./(Gas_density(1)*M*cp_G); %Coeff. dispersione termica m^2/s
        alfa(i) = E_t(i)/Interface_gasDensity/cp_G;
        delta_Speed = Speed_array(i) - (v_star+alfa(i)./Temperature_array(i).*(dTdy(i)-dTdy(end)));
        if delta_Speed*delta_Speed_LL < 0
            Speed_UL = Speed_array(i);
            
        else
            Speed_LL = Speed_array(i);
            delta_Speed_LL = delta_Speed;
        end
    end
end
dTdt = zeros((NumbersSpatialDiscretisation),1);
dCdt = dTdt;
for i = 2 : (NumbersSpatialDiscretisation-1)
    % calculations of the time derivative
    dTdt_i = (v_avg .* y(i) + (v_star - v_avg)) ./ (v_avg .* t - H) .* dTdy(i) + (alfa(i) ./ (Temperature_array(1) .* (v_avg .* t - H).^2)) * (Temperature_array(i).* d2Tdy2(i) - dTdy(i).^2 + dTdy(1) .* dTdy(i));
    dCdt_i =(v_avg .* y(i) + (v_star - v_avg)) ./ (v_avg .* t - H) .* dCdy(i) + (E_m(i)./ (v_avg .* t - H).^2) .* d2Cdy2(i) - alfa(i)./ (Temperature_array(1).*(v_avg .* t - H).^2) .* (dCdy(i).* (dTdy(i) - dTdy(1)) + Concentration_array(i).* d2Tdy2(i));
    dTdt(i) = dTdt_i;
    dCdt(i) = dCdt_i;
    %  % corrections if something impossibile for a physical POV happens
    % 
    %  % the conditions for dTdt(i) being zeroed are having a point
    %  % temperature higher than the roof temperature and it not decreasing
    %  % i.e. its time partial derivative is positive
    % 
    %  if (Temperature_array(i)) > Temperature_array(end) && dTdt_i>0
    %      %dTdt(i) = (v_avg .* y(i) + (v_star - v_avg)) ./ (v_avg .* t - H) .* dTdy(i) + (alfa(i) ./ (Temperature_array(1) .* (v_avg .* t - H).^2)) * (Temperature_array(i).* d2Tdy2(i) - dTdy(i).^2 + dTdy(1) .* dTdy(i));
    %      dTdt(i) = 0;
    % 
    %  else
    %      dTdt(i) = dTdt_i;
    % 
    % end
    % 
    %  % the conditions for dCdt(i) being zeroed are having a point
    %  % concentration lower than 0 and it not increasing
    %  % i.e. its time partial derivative is negative
    % 
    %  if (Concentration_array(i)) < 0 && dCdt_i<0
    % 
    % 
    %      dCdt(i) = 0;
    %  else
    %      dCdt(i) = dCdt_i;
    % 
    % 
    %  end
    %dCdt(i) =(v_avg .* y(i) + (v_star - v_avg)) ./ (v_avg .* t - H) .* dCdy(i) + (E_m(i)./ (v_avg .* t - H).^2) .* d2Cdy2(i) - alfa(i)./ (Temperature_array(1).*(v_avg .* t - H).^2) .* (dCdy(i).* (dTdy(i) - dTdy(1)) + Concentration_array(i).* d2Tdy2(i));
end
%Dirichelet e Neumann in y=1
% Temperature_array(end)=T_t;
if Speed_array(end)>=0
    % BC for filling and rest: zero space partial derivative at roof
    % spatial coordinate
    Concentration_array_end=Concentration_array(end-1);
    Temperature_array_end = T_t;
else
    % BC for emptying: flow conservation for concentration, IMO should be
    % the case for energy too, but does not work. Needs to be recalculated
    % not in 15 ms, needs some thinking. Also assuming that T = T_inletGas
    % could work. Something [redacted] off whenever emptying occurs, in point
    % directly before the roof something is off, needs to be looked at.
    Concentration_array_end = Concentration_array(end-1) + Concentration_array(end)*Speed_array(end)*(H-v_avg*t)/E_m(end)*dy;
    %Temperature_array_end = Temperature_array(end-1) + Temperature_array(end)/cp_G*Speed_array(end)*(H-v_avg*t)/E_t(end)*dy;
    Temperature_array_end = T_t;
end
Output_array = zeros(2*NumbersSpatialDiscretisation,1);
% defining and filling the output array for the ODE
Output_array(1) = Concentration_array(1) - SaturationConcentration;
for i = 2 : (NumbersSpatialDiscretisation-1)
    Output_array(i) = dCdt(i);
end
Output_array(NumbersSpatialDiscretisation) = Concentration_array_end - Concentration_array(NumbersSpatialDiscretisation);
%Output_array(NumbersSpatialDiscretisation+1) = calc_T_ILG(Temperature_array(1),TL, dTdy(1), dCdy(1), v_avg, h_L, k_G, cpL, rho_liq, DeltaH_evap);
Output_array(NumbersSpatialDiscretisation+1) = TL - Temperature_array(1);
for i = 2 : (NumbersSpatialDiscretisation-1)
    Output_array(NumbersSpatialDiscretisation + i) = dTdt(i);
end
Output_array(2*NumbersSpatialDiscretisation) = Temperature_array_end - Temperature_array(end);
end
function phi = flux_limiter(r, method)
switch method
    case 'minmod'
        phi = max(0, min(1, r));  % Minmod limiter
    case 'superbee'
        phi = max(0,max(min(2*r, 1), min(r, 2))); % Superbee limiter
    case 'vanleer'
        phi = (r + abs(r)) / (1 + abs(r));  % Van Leer limiter
    otherwise
        phi = 1;  % No limiter (upwind scheme)
end
end
%% Funzione per calcolare v_avg durante i cicli
function v_avg = v_cicli(t, V_serbatoio, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration)
v_riempimento = V_serbatoio / (pi * (diametro / 2)^2 * t_riempimento);
v_svuotamento = -V_serbatoio / (pi * (diametro / 2)^2 * t_svuotamento);
cycle_time = mod(t, cycle_duration);
if cycle_time <= t_riempimento
    v_avg = v_riempimento;
elseif cycle_time <= (t_riempimento + t_riposo)
    v_avg = 0;
else
    v_avg = v_svuotamento;
end
end

This is the new one:

%% New code
clear;
clc;
global  M R_gas k D cp_G mu_G cpL h_L rho_liq DeltaH_evap H_c TL A_ant B_ant C_ant NumbersSpatialDiscretisation t_riempimento t_riposo t_svuotamento cycle_duration
%% inputs 
T_t = 328;  %268,283,288,298,328
%% Parametri del serbatoio e costanti fisiche
H = 15; % Altezza del serbatoio (m)
diametro = H; % Diametro del serbatoio (m)
V_serbatoio = pi * (diametro / 2)^2 * (H-1); % Volume del serbatoio (m^3)
P = 1e5; % Pressione in Pa
%% Costanti Diesel
M = 0.14; % Massa molare del Diesel in kg/mol
R_gas = 8.314; % Costante dei gas perfetti in J/(mol*K)
k = 0.025; % Conducibilità termica dell'aria in W/m·K
D = 10^(-5); % Diffusività in m²/s
cp_G = 2100; % Calore specifico in J/kg·K
mu_G = 1.85e-5; % Viscosità dinamica dell'aria (Pa·s)
cpL = 2200; % Calore specifico Diesel liquido in J/kg·K
h_L = 12; % Coefficiente scambio termico del liquido W/m²·K
rho_liq = 835; % Densità del Diesel liquido a 298 K in kg/m³
DeltaH_evap = 375000*M; % Calore latente di evaporazione in J/mol
H_c = 0; % Costante di Henry NON UTILIZZATA
TL = 283; % Temperatura liquido in K
%% Costanti Antoine (per tensione di vapore)
A_ant = 12.101;
B_ant = 8907;
C_ant =0;
%% Parametri discretizzazione
NumbersSpatialDiscretisation = 50; % Numero di nodi spaziali
%% Condizioni iniziali
T_init = T_t * ones(NumbersSpatialDiscretisation, 1); % Temperatura iniziale uniforme (K)
C_init = 0.01*ones(NumbersSpatialDiscretisation,1);
%C_init = linspace(0.02, 0.01, NumbersSpatialDiscretisation)'; % Gradiente iniziale della concentrazione
y0 = [C_init; T_init]; % Condizioni iniziali combinate
%% Parametri per i cicli di riempimento, riposo e svuotamento
t_riempimento = 24 * 3600; % 24 ore
t_riposo = 1 * 24 * 3600; % 5 giorni di riposo
t_svuotamento = 24 * 3600; % 24 ore di svuotamento
cycle_duration = t_riempimento + t_riposo + t_svuotamento;
transition_time = 3600; % Tempo di transizione (1 ora)
%% Funzione per calcolare v_avg
v_avg_func = @(t) v_cicli(t, V_serbatoio, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time);
%% Definizione dei parametri extra da passare
params = {NumbersSpatialDiscretisation, D, k, cp_G, mu_G, diametro, v_avg_func, P, M, R_gas, h_L, rho_liq, cpL, DeltaH_evap, A_ant, B_ant, C_ant, H, H_c,V_serbatoio, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time,TL,T_t};
%% Matrice di massa (per condizioni al contorno)
mass_matrix = eye(2 * NumbersSpatialDiscretisation);
mass_matrix(1, 1) = 0;
mass_matrix(NumbersSpatialDiscretisation, NumbersSpatialDiscretisation) = 0;
mass_matrix(NumbersSpatialDiscretisation + 1, NumbersSpatialDiscretisation + 1) = 0;
mass_matrix(2 * NumbersSpatialDiscretisation, 2 * NumbersSpatialDiscretisation) = 0;
t_span = [0,cycle_duration*30]
%% Risoluzione del sistema tramite ode15s per ogni intervallo temporale
Output = OdeCompiler(t_span, params{:},mass_matrix,y0);
t = Output(:,1);
u_sol = Output(2:end,:);
%% FCN with time related jumps
function out = OdeCompiler(t_span, ~, ~, ~, ~, ~, diametro, ~, P, ~, ~, ~, ~, ~, ~, ~, ~, ~, H, ~,V_serbatoio, ~, ~, ~, ~, ~,~,T_t,~,~)
global  M R_gas k D cp_G mu_G cpL h_L rho_liq DeltaH_evap H_c TL A_ant B_ant C_ant NumbersSpatialDiscretisation t_riempimento t_riposo t_svuotamento cycle_duration
%% Parametri discretizzazione
y = linspace(0,1,NumbersSpatialDiscretisation);
%% Condizioni iniziali
T_init = TL * ones(NumbersSpatialDiscretisation, 1); % Temperatura iniziale uniforme (K)
T_init(NumbersSpatialDiscretisation) = T_t;
C_init = 0.05*ones(NumbersSpatialDiscretisation,1);
%C_init = linspace(0.02, 0.01, NumbersSpatialDiscretisation)'; % Gradiente iniziale della concentrazione
y0 = [C_init; T_init]; % Condizioni iniziali combinate
%% Parametri per i cicli di riempimento, riposo e svuotamento
t_riempimento = 12 * 3600; % 12 ore
t_riposo = 10 * 24 * 3600; % 7 giorni di riposo
t_svuotamento = 20*24 * 3600; % 24 ore di svuotamento
cycle_duration = t_riempimento + t_riposo + t_svuotamento+t_riposo/30;
transition_time = 3600; % Tempo di transizione (1 ora)
%% Funzione per calcolare v_avg
v_avg_func = @(t) v_cicli(t, V_serbatoio, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time);
%% Definizione dei parametri extra da passare
params = {NumbersSpatialDiscretisation, D, k, cp_G, mu_G, diametro, v_avg_func, P, M, R_gas, h_L, rho_liq, cpL, DeltaH_evap, A_ant, B_ant, C_ant, H, H_c,V_serbatoio, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time,TL,T_t};
%% Matrice di massa (per condizioni al contorno)
mass_matrix = eye(2 * NumbersSpatialDiscretisation);
mass_matrix(1, 1) = 0;
mass_matrix(NumbersSpatialDiscretisation, NumbersSpatialDiscretisation) = 0;
mass_matrix(NumbersSpatialDiscretisation + 1, NumbersSpatialDiscretisation + 1) = 0;
mass_matrix(2 * NumbersSpatialDiscretisation, 2 * NumbersSpatialDiscretisation) = 0;
%compilazione del vettore tempo dividendolo in range con velocità del
%fluido costante per evitare cambi improvvisi durente la risoluzione ODE
NumberOfCycles = ceil(t_span(2)/cycle_duration);
Time_vector = zeros(NumberOfCycles*3+1,1);
for i = 0: (NumberOfCycles-1)
    Time_vector(i*3+2) = cycle_duration*i + t_riempimento;
    Time_vector(i*3+3) = cycle_duration*i + t_riempimento + t_riposo;
    Time_vector(i*3+4) = cycle_duration*i + t_riempimento + t_riposo +  t_svuotamento;
end
v_riempimento = V_serbatoio / (pi * (diametro / 2)^2 * t_riempimento);
v_svuotamento = -0.3*V_serbatoio / (pi * (diametro / 2)^2 * t_svuotamento);
timeSolution_vector = zeros(1,1);
Solution_vector = zeros(2,NumbersSpatialDiscretisation);
% calcolo velocità del fluido
for i = 2 : length(Time_vector)
    cycle_time = mod(Time_vector(i), cycle_duration+0.001);
    if cycle_time <= t_riempimento
        v_avg = v_riempimento;
    elseif cycle_time <= (t_riempimento + t_riposo)
        v_avg = 0;
    else
        v_avg = v_svuotamento;
    end
    % parametri ODE
    params = {NumbersSpatialDiscretisation, D, k, cp_G, mu_G, diametro, v_avg_func, P, M, R_gas, h_L, rho_liq, cpL, DeltaH_evap, A_ant, B_ant, C_ant, H, H_c,V_serbatoio, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time,TL,T_t,v_avg};
    tic 
    % Define options for ode15s
    opts = odeset('RelTol', 1e-10, 'AbsTol', 1e-10, 'Mass', mass_matrix, 'MassSingular', 'yes');
    
    if i > 2
        disp(i)
        pde_fun(t(end), u_sol(end,:)', params{:});
    else
    end
    % Solve the DAE system
    [t, u_sol] = ode15s(@(t, u) pde_fun(t, u, params{:}), t_span, y0, opts);
    % Post-process the solution if needed
    toc
    if i == 2
        Solution_vector = u_sol;
        timeSolution_vector = t;
    else
        Solution_vector = [Solution_vector; u_sol];
        timeSolution_vector = [timeSolution_vector;t];
    end
    y0 = u_sol(end,:)';
    
    pde_fun(t(1), u_sol(1,:)', params{:});
end
out = [timeSolution_vector,Solution_vector];
figure(1)
plot(timeSolution_vector/3600,Solution_vector(:,NumbersSpatialDiscretisation))
title('Andamento Ctop nel tempo')
xlabel('tempo [h]')
ylabel('Ctop [mol/m3]')
figure(2)
plot(y,Solution_vector(2000,1:NumbersSpatialDiscretisation))
title('Andamento C nello spazio')
xlabel('y[-]')
ylabel('C [mol/m3]')
figure(3)
plot(y,Solution_vector(2000,NumbersSpatialDiscretisation+1:2*NumbersSpatialDiscretisation))
title('Andamento T nello spazio')
xlabel('y[-]')
ylabel('T [K]')
figure(4)
surf(y,timeSolution_vector/3600,Solution_vector(:,1:NumbersSpatialDiscretisation))
title('grafico 3D C')
xlabel('y[-]')
ylabel('tempo [h]')
zlabel('C [mol/m3]')
figure(5)
surf(y,timeSolution_vector/3600,Solution_vector(:,NumbersSpatialDiscretisation+1:2*NumbersSpatialDiscretisation))
title('grafico 3D T')
xlabel('y[-]')
ylabel('tempo [h]')
zlabel('T [K]')
pde_fun(timeSolution_vector(1), Solution_vector(1,:),params{:});
end
%%funzione PDE
function Output_array = pde_fun(t, u, NumbersSpatialDiscretisation, D, k_G, cp_G, mu_G, diametro, ~, P, M, R_gas, ~, ~, ~, ~, A_ant, B_ant, ~, H, ~,~, t_riempimento, t_riposo, ~, cycle_duration, ~,TL,T_t,v_avg)
% Adaptation in time for solver stability during filling and emptying
number_cycles = floor(t / cycle_duration);
if v_avg < 0
    t = t - (t_riempimento + t_riposo) * number_cycles;
elseif v_avg > 0
    t = t - cycle_duration * number_cycles;
end
% Receiving input
Concentration_array = u(1:NumbersSpatialDiscretisation);
Temperature_array = u((NumbersSpatialDiscretisation + 1):(2 * NumbersSpatialDiscretisation));
% Spatial derivative discretization
dTdy = zeros(NumbersSpatialDiscretisation, 1);
d2Tdy2 = zeros(NumbersSpatialDiscretisation, 1);
dCdy = zeros(NumbersSpatialDiscretisation, 1);
d2Cdy2 = zeros(NumbersSpatialDiscretisation, 1);
dy = 1 / NumbersSpatialDiscretisation;
y = linspace(0, 1, NumbersSpatialDiscretisation);
for i = 2:(NumbersSpatialDiscretisation - 1)
    % Concentration derivative with flux limiter
    rC = (Concentration_array(i) - Concentration_array(i - 1)) / (Concentration_array(i + 1) - Concentration_array(i));
    limiterC = flux_limiter(rC, 'minmod');
    dCdy(i) = limiterC * (Concentration_array(i + 1) - Concentration_array(i)) / dy;
    % Temperature derivative with flux limiter
    rT = (Temperature_array(i) - Temperature_array(i - 1)) / (Temperature_array(i + 1) - Temperature_array(i));
    limiterT = flux_limiter(rT, 'minmod');
    dTdy(i) = limiterT * (Temperature_array(i + 1) - Temperature_array(i)) / dy;
    % Second derivatives using central difference
    d2Tdy2(i) = (Temperature_array(i + 1) - 2 * Temperature_array(i) + Temperature_array(i - 1)) / dy^2;
    d2Cdy2(i) = (Concentration_array(i + 1) - 2 * Concentration_array(i) + Concentration_array(i - 1)) / dy^2;
end
dCdy(1) = (Concentration_array(2) - Concentration_array(1)) / dy;
dTdy(NumbersSpatialDiscretisation) = (Temperature_array(NumbersSpatialDiscretisation) - Temperature_array(NumbersSpatialDiscretisation - 1)) / dy;
dCdy(NumbersSpatialDiscretisation) = (Concentration_array(NumbersSpatialDiscretisation) - Concentration_array(NumbersSpatialDiscretisation - 1)) / dy;
% Raoult's law for saturation concentration assuming ideal gas
SaturationConcentration = 6894.76 * exp(A_ant - B_ant / (1.8 * Temperature_array(1))) / P * P / (R_gas * Temperature_array(1));
Interface_gasDensity = P / (R_gas * Temperature_array(1));
% Define v*
v_star = v_avg - (D / Interface_gasDensity) * dCdy(1);
% Initialize iterative calculation for transport coefficients
Speed_array = zeros(NumbersSpatialDiscretisation, 1);
Gas_density = Speed_array;
Re = Speed_array;
Sc = Speed_array;
Pe_m = Speed_array;
Pe_t = Speed_array;
E_m = Speed_array;
E_t = Speed_array;
alfa = Speed_array;
Gas_density(1) = Interface_gasDensity;
SpeedPreviousIteration = v_star;
LowerLimit = SpeedPreviousIteration / 5;
% Upper limit calculations
UpperLimit = SpeedPreviousIteration * 5;
Re(1) = UpperLimit * diametro * Gas_density(1) / mu_G;
Re(1) = abs(Re(1));
Pr = mu_G * cp_G / k_G;
Sc(1) = mu_G / (Gas_density(1) * M * D);
Pe_m(1) = (1 / (Re(1) * Sc(1)) + (Re(1) * Sc(1)) / 192);
Pe_t(1) = (1 / (Re(1) * Pr) + (Re(1) * Pr) / 192);
E_m(1) = Pe_m(1) * D;
E_t(1) = Pe_t(1) * k_G / (Gas_density(1) * M * cp_G);
alfa(1) = E_t(1) / Interface_gasDensity / cp_G;
% Initialize delta_Speed_LL
delta_Speed_LL = -inf;
% Bisection method for true vapor speed at the interface
Speed_array(1) = LowerLimit + (UpperLimit - LowerLimit) / 2;
delta_Speed = Speed_array(1) - (v_star + alfa(1) / Temperature_array(1) * (dTdy(1) - dTdy(end)));
itCount = 0;
while abs(delta_Speed) / abs(Speed_array(1)) > 1e-3 && itCount < 100
    itCount = itCount + 1;
    Speed_array(1) = LowerLimit + (UpperLimit - LowerLimit) / 2;
    Re(1) = Speed_array(1) * diametro * Gas_density(1) / mu_G;
    Re(1) = abs(Re(1));
    Sc(1) = mu_G / (Gas_density(1) * M * D);
    Pe_m(1) = (1 / (Re(1) * Sc(1)) + (Re(1) * Sc(1)) / 192);
    Pe_t(1) = (1 / (Re(1) * Pr) + (Re(1) * Pr) / 192);
    E_m(1) = Pe_m(1) * D;
    E_t(1) = Pe_t(1) * k_G / (Gas_density(1) * M * cp_G);
    alfa(1) = E_t(1) / Interface_gasDensity / cp_G;
    delta_Speed = Speed_array(1) - (v_star + alfa(1) / Temperature_array(1) * (dTdy(1) - dTdy(end)));
    if delta_Speed * delta_Speed_LL < 0
        UpperLimit = Speed_array(1);
    else
        LowerLimit = Speed_array(1);
        delta_Speed_LL = delta_Speed;
    end
end
for i = 2:NumbersSpatialDiscretisation
    % Upper limit calculations
    Speed_UL = Speed_array(1) * 5;
    Gas_density(i) = P / (R_gas * (TL + T_t) / 2);
    Re(i) = Speed_UL * diametro * Gas_density(i) / mu_G;
    Re(i) = abs(Re(i));
    Sc(i) = mu_G / (Gas_density(i) * M * D);
    Pe_m(i) = (1 / (Re(i) * Sc(i)) + (Re(i) * Sc(i)) / 192);
    Pe_t(i) = (1 / (Re(i) * Pr) + (Re(i) * Pr) / 192);
    E_m(i) = Pe_m(i) * D;
    E_t(i) = Pe_t(i) * k_G / (Gas_density(1) * M * cp_G);
    alfa(i) = E_t(i) / Interface_gasDensity / cp_G;
    % Lower limit calculations
    Speed_LL = Speed_array(1) / 5;
    Gas_density(i) = P / (R_gas * Temperature_array(i));
    Re(i) = Speed_LL * diametro * Gas_density(i) / mu_G;
    Re(i) = abs(Re(i));
    Sc(i) = mu_G / (Gas_density(i) * M * D);
    Pe_m(i) = (1 / (Re(i) * Sc(i)) + (Re(i) * Sc(i)) / 192);
    Pe_t(i) = (1 / (Re(i) * Pr) + (Re(i) * Pr) / 192);
    E_m(i) = Pe_m(i) * D;
    E_t(i) = Pe_t(i) * k_G / (Gas_density(1) * M * cp_G);
    alfa(i) = E_t(i) / Interface_gasDensity / cp_G;
    delta_Speed_LL = Speed_LL - (v_star + alfa(i) / Temperature_array(i) * (dTdy(i) - dTdy(end)));
    % Bisection method for true speed at a given point
    Speed_array(i) = Speed_LL + (Speed_UL - Speed_LL) / 2;
    delta_Speed = Speed_array(i) - (v_star + alfa(i) / Temperature_array(i) * (dTdy(i) - dTdy(end)));
    itCount = 0;
    while abs(delta_Speed / Speed_array(1)) > 1e-3 && itCount < 100
        itCount = itCount + 1;
        Speed_array(i) = Speed_LL + (Speed_UL - Speed_LL) / 2;
        delta_Speed = Speed_array(i) - (v_star + alfa(i) / Temperature_array(i) * (dTdy(i) - dTdy(end)));
        if delta_Speed * delta_Speed_LL < 0
            Speed_UL = Speed_array(i);
        else
            Speed_LL = Speed_array(i);
            delta_Speed_LL = delta_Speed;
        end
    end
end
dTdt = zeros(NumbersSpatialDiscretisation, 1);
dCdt = dTdt;
for i = 2:(NumbersSpatialDiscretisation - 1)
    dTdt_i = (v_avg * y(i) + (v_star - v_avg)) / (v_avg * t - H) * dTdy(i) + (alfa(i) / (Temperature_array(1) * (v_avg * t - H)^2)) * (Temperature_array(i) * d2Tdy2(i) - dTdy(i)^2 + dTdy(1) * dTdy(i));
    dCdt_i = (v_avg * y(i) + (v_star - v_avg)) / (v_avg * t - H) * dCdy(i) + (E_m(i) / (v_avg * t - H)^2) * d2Cdy2(i) - alfa(i) / (Temperature_array(1) * (v_avg * t - H)^2) * (dCdy(i) * (dTdy(i) - dTdy(1)) + Concentration_array(i) * d2Tdy2(i));
    dTdt(i) = dTdt_i;
    dCdt(i) = dCdt_i;
    if Temperature_array(i) > Temperature_array(end) && dTdt_i > 0
        dTdt(i) = 0;
    end
    if Concentration_array(i) < 0 && dCdt_i < 0
        dCdt(i) = 0;
    end
end
if Speed_array(end) >= 0
    Concentration_array_end = Concentration_array(end - 1);
    Temperature_array_end = T_t;
else
    Concentration_array_end = Concentration_array(end - 1) + Concentration_array(end) * Speed_array(end) * (H - v_avg * t) / E_m(end) * dy;
    Temperature_array_end = T_t;
end
Output_array = zeros(2 * NumbersSpatialDiscretisation, 1);
Output_array(1) = Concentration_array(1) - SaturationConcentration;
for i = 2:(NumbersSpatialDiscretisation - 1)
    Output_array(i) = dCdt(i);
end
Output_array(NumbersSpatialDiscretisation) = Concentration_array_end - Concentration_array(NumbersSpatialDiscretisation);
Output_array(NumbersSpatialDiscretisation + 1) = TL - Temperature_array(1);
for i = 2:(NumbersSpatialDiscretisation - 1)
    Output_array(NumbersSpatialDiscretisation + i) = dTdt(i);
end
Output_array(2 * NumbersSpatialDiscretisation) = Temperature_array_end - Temperature_array(end);
end
function phi = flux_limiter(r, method)
switch method
    case 'minmod'
        phi = max(0, min(1, r));
    case 'superbee'
        phi = max(0, max(min(2 * r, 1), min(r, 2)));
    case 'vanleer'
        phi = (r + abs(r)) / (1 + abs(r));
    otherwise
        phi = 1;
end
end
%% Funzione per calcolare v_avg durante i cicli
function v_avg = v_cicli(t, V_serbatoio, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration)
v_riempimento = V_serbatoio / (pi * (diametro / 2)^2 * t_riempimento);
v_svuotamento = -V_serbatoio / (pi * (diametro / 2)^2 * t_svuotamento);
cycle_time = mod(t, cycle_duration);
if cycle_time <= t_riempimento
    v_avg = v_riempimento;
elseif cycle_time <= (t_riempimento + t_riposo)
    v_avg = 0;
else
    v_avg = v_svuotamento;
end
end

[Walter Roberson]

To be honest, your code has a bad case of "tl;dr" (Too Long, Didn't Read). It would take notably more work to analyze then I am willing to put in.

[Diego Rampi]

You're right, the code is very complex, partly because it addresses a complex problem itself, and partly because, being self-taught, I’m not exactly a champion of optimization and readability. I’ll try to explain the problem conceptually:

The issue is that while I previously had these constant values, they are now provided to me from an Excel file as discrete variables. Ode15s, to solve the problem, wants to adapt the discretization as it sees fit, and when it needs a value at a time point not present in the Excel file, it gets stuck. There are two possible solutions:

1. Transform the two discrete inputs into continuous functions (this seems very difficult because they’re really messy). 2. Tell ode15s that when it doesn’t find the exact value, it can use the closest one to the required time point (or the result of an interpolation).

I hope I was clear. Thank you!

Answer 1

[Walter Roberson]

interp1() with default interpolation, and interp1() 'nearest' seldom have a place in ode*() calls. The problem is that the interpolated values hae discontinuous first derivatives at every breakpoint, but the mathematics of ode*() requires that all values have continuous second derivative. interp1('spline') generally satisfies the mathematics required.

Whether interp1('spline') satisfies the underlying problem is a different question. If the values being interpolated are understood as being discrete samplings of continuous functions, then spline interpolation is generally fine. But if the values being sampled are understood as representing discrete events, then spline interpolation is not suitable.

If for example the data represents times that chemicals are injected into a system, then the chemicals are absent over time t1 to t2 and then suddenly appear at t2 sharp. Even if they are injected over time starting at t2, it is the case that they are absent from t1 to t2 and then change -- a flat derivative from t1 to t2 and then suddenly a different derivative at t2. A spline interpolation of values from t1 to t2 would misrepresent the system.

If you have discrete events happening, then you need to stop the ode*() at each breakpoint, adjust the boundary conditions, and restart the ode*() for the remaining time. If the discrete events happen at precise times then the easiest way to handle this situation is to specify ode*() tspan that end at the next discrete event; if the discrete events happen at internal boundaries (e.g., ball bouncing against a wall) then you need to use ode event functions.