I’m trying to solve for a vector k (reaction constant) that gets updated every time T is updated, and that new value of k is used to find the next T and so on. But the first time T is called for each iteration of k should use the value of T corresponding to the same index of k, and the second time T is called (i.e. T_prev) should reference the value of the previous index. I’m not really sure what I’m doing wrong. I usually don’t use function handles unless I have to.

NA_0 = 3170; % mols of ONCB
NB_0 = 43000; % mols of NH3
B = NB_0 / NA_0;
dHrx = -5.9E5;
Vn = 3.265; % Volume of reactant ONCB at normal operating conditions
T_a = 298; % ambient temperature (K)

C_pA = 40E-3; % heat capacity of ONCB
C_pB = 8.38E-3; % heat capacity of H2O
C_pW = 18E-3; % heat capacity of NH3

UA = 35.85;
xt = linspace(0, 120, 240)'; % time vector

% Initial rate constant k at T_a
k = 1.167E-4;

tspan = [0 120]; 
ic = 448; % initial temperature in K

% ODE function for temperature change
[t_s, T] = ode15s(@(t_s, T) myode(t_s, T, NA_0oc, NB_0oc, B, dHrx, Vn, T_a, C_pA, C_pB, C_pW, UA), tspan, ic);

% Plot the results
plot(t_s, T)
title("Excess Reagents")
xlabel('Time (minutes)')
ylabel('Temperature (K)')
title('Temperature vs Time')


function dTdt = myode(t, T, NA_0, NB_0, B, dHrx, Vn, T_a, C_pA, C_pB, C_pW, UA)
    % Calculate the current rate constant k based on temperature
    
    if T == ic
        T_prev = ic - 0.1;
        k = @(T, T_prev) 1.167E-4 * exp(11273 / (1.987 * ((1/T) - (1 ./ T_prev))));
    else
        for z = 2:length(T)
            T_prev(z) = T(z-1);
            k = @(T, T_prev) 1.167E-4 * exp(11273 / (1.987 * ((1/T) - (1 ./ T_prev))));
        end
    end
  
    % Calculate dX/dt
    xt = NA_0 - NB_0 * (1 - exp(-k(T, T_prev) * t)); % Assuming first-order kinetics
    dXdt = k(T, T_prev) * NA_0 .* (1 - xt) .* (B - 2 .* xt) ./ Vn;
    
    % Calculate the heat generated by the reaction
    Q_b = k(T, T_prev) * (NA_0^2) .* (1 - xt) .* (B - 2 .* xt) .* (-dHrx) / Vn;

    % Calculate the heat removed by the reactor cooling
    Q_r = UA * (T - T_a);

    % Calculate the total heat capacity
    NC_p = NA_0 * C_pA + NB_0 * C_pB + NB_0 * C_pW;

    % Check if any element of Q_b is infinite
    if any(isinf(Q_b))
        dTdt = 0;
    else
        % Calculate the rate of temperature change
        dTdt = (Q_b + Q_r) / NC_p; % Heat generated - Heat removed
    end
end

This is the error I get currently:

Unrecognized function or variable ‘k’.

Error in Monsanto Plant Disaster>myode (line 61)
xt = NA_0 – NB_0 * (1 – exp(-k(T, T_prev) * t)); % Assuming first-order kinetics

Error in Monsanto Plant Disaster (line 24) [t_s, T] = ode15s(@(t_s, T)
myode(t_s, T, NA_0, NB_0, B, dHrx, Vn, T_a, C_pA, C_pB, C_pW, UA),
tspan, ic);

Error in odenumjac (line 131)
Fdel(:,j) = F(Fargs{1:diffvar-1},ydel(:,j),Fargs{diffvar+1:end});

Error in ode15s (line 349)
[dfdy,Joptions.fac,nF] = odenumjac(ode, {t,y,odeArgs{:}}, f0, Joptions); %#ok