Why is my state.u matrix a constant and how to make my heat time dependent?

Basic science behind my code:

I have my internal heat generation dependent on two parameters - alpha (degree of hydration) and temperature. Heat generation is a function of dalpha/dt which is inturn a function of current alpha and temperature. alpha is updated for every iteration by multiplying dalpha/dt with time step. So, in short, alpha is just a function of temperature and time and so my whole heat generation is a function of temperature and time.

My equations:

=

where all are constants except kd, De, kr which are temperature dependent.

So for every iteration, we calculate alpha from previous dalpha/dt use it in the heat equation, calculate the new temperature, then alpha, then dalpha/dt and so on.

That means my heat generation is dependent on both temperature and time.

My question is:

How to make this code take the instantaneous temperature in calculating heat generation? So far it is working for a constant initial temperature. When I have gone through state.u matrix, it is just a row of initial temperature conditions..so how to update state.u conditions and use the updated state.u values in heat generation iterations?

function q = myfun(~, state) 
if isnan(state.u)
    q = NaN(size(state.u));
    return;
end
if isnan(state.time)
    q = NaN(size(state.u));
    return;
end
else
%Parameters
Rho_ce = 3.15; %g/cm3 == density of cement
Rho_w = 1; %g/cm3 == density of water
v_cbw = 0.25; %gmass of chemically bound water/gcement
wg = 0.15; %gmass of physically bound water/gcement
r0 = 9*10^-4;%cm
T_initial = 303;
wc_ratio = 0.303;% Water to cement ration
k_3 = (0.4/wc_ratio);

C3S = 54;C2S = 19;
C3A = 11;C4AF = 10;
Gypsum = 6;

beta_1 = 1000;beta_2 = 1000;
E_R = 5400;beta_3 = 7500;

Cem = 565*10^-9;
He = 5.9444*10^5;
tstep = 100;


M = length(state.u);
alpha = real(zeros(1, M));
dalphadt = real(zeros(1, M));

R3 = zeros(1, M);S3 = zeros(1, M);
T3 = zeros(1, M);U3 = zeros(1, M);
q = real(zeros(1, M));

R3(:,1) = (((6*10^-12)*(C3S + C3A)) + (4*10^-10))*exp(-beta_1*((1/T_initial)-(1/293)));
S3(:,1) = (0.0003*C3S + 0.0186)*exp(-beta_2*((1/T_initial)-(1/293)));
T3(:,1) = r0./((7*10^-10 - (8*10^-12)*C2S)*(exp(-beta_3*((1/T_initial)-(1/293)))));
U3(:,1) = 1./(((8*10^-8)*C3S + (10^-6))*exp(-E_R*((1/T_initial)-(1/293))));

for i = 2:M
        alpha(1, i) = alpha(1, i-1) + dalphadt(1, i-1).*(tstep/3600);
        R3(1,i) = (((6*10^-12)*(C3S + C3A)) + (4*10^-10))*exp(-beta_1.*((1./state.u(1, i))-(1/293)));
        S3(1,i) = (0.0003*C3S + 0.0186)*exp(-beta_2.*((1./state.u(1, i))-(1/293)));
        T3(1,i) = r0./((7*10^-10 - (8*10^-12)*C2S)*(exp(-beta_3.*((1./state.u(1, i))-(1/293)))));
        U3(1,i) = 1./(((8*10^-8)*C3S + (10^-6))*exp(-E_R.*((1./state.u(1, i))-(1/293))));
        dalphadt(1, i) = ((3*Rho_w/(r0*Rho_ce*(v_cbw+wg))) .* power(1 - k_3 .* alpha(1, i), (2.6-4*wc_ratio))) ./ ((1 ./ ((R3(1, i) .* power(alpha(1, i), -1.5)) + S3(1, i) .* power(alpha(1, i), 3))) + T3(1, i) ./ log(alpha(1, i)) + U3(1, i) .* power(1 - alpha(1, i), -2/3) - T3(1, i) .* power(1 - alpha(1, i), -1/3) ./ log(alpha(1, i)));    
        q(1, i) = (Cem*He).*dalphadt(1, i)./3600;
        
end
end
end

[Kevin Isakayoga]

Up! I also want to know.