Numerical solution of an iterative equation

I am trying to solve the following iterative equation:

where the parameter tao (needed in the previous equation) is defined as:

Here’s the parameters and my approach to solving the equation in Mathematica:

ah = 342496; (*J/mol*)
R = 8.314; (*J/mol.K*)
A = 7.6*10^-38; (*s*)
b = 0.67;
x = 0.49;
q = -0.1/60;(*K/s*)
T0 = 500;(*K*)
Tfinal = 350;(*K*)
dt = 0.1 ;(*s*)
n = (T0 - Tfinal)/dt; (*number of steps*)
Tf0 = T0;
dT = dt*q; (*K*)

Tk = T0 + Sum[dT, {i, 1, n}]; 
τ = A*Exp[(x ah)/(R Tk) + (1 - x) (ah/Subscript[T, fkminus1])]; 
Tf = T0 + Sum[dT (1 - Exp[-Sum[(dT/(q*τ))^b, {i, 1, n}]]), {i, 1, n}]

In the code, notice that dT = dt*q, where dt is a time step that I choose which ideally should be small enough as to ensure that Tf decays less than 1 K. Also, notice that $tau_{0,k}$ is a function of $T_{f,k-1}$ ($T_{f,0}$ is defined in the code as 500) and $T_k$.

The problem I am having is that I do not know how to connect the two equations so that the $T_k$ that I calculate from equation 1, eventully goes to equation 2 and the $tau_k$ that I calculate from equation 2 goes to equation 1.

Questions

1. How can I solve the Tf equation?
2. How can I plot Tf vs T, similar to the figure below (only for one q, in this case for q = -0.1/60)?. Notice that the parameters from figure below are the same I am using in the code.

Note: Here’s the paper where the two equations are coming from as well as the figure (figure 1 in that paper) in case anybody wants to look at it: https://sci-hub.se/10.1016/s0260-8774(02)00212-1.

I hope everything is clear and I appreciate in advance your help.

EDIT: Clarification on the parameters

1. ΔT is the change in temperature and it is related by the change in time (which is set by the user, small enough so that Tf doesn’t decrease more than 1 K) as ΔT=dt*q. Both, dt and q can be taken as constant. ΔTk represents then a series of temperature steps.
2. qk is the cooling rate in Kelvin/seconds at which the process is performed. In the code it fixed with value of -0.1/60 (or -0.0016 Kelvin/second). The negative sign is because it is a cooling process (it would be positive if it was heating).
3. T,k is simply the temperature, which is defined as Tk = T0 + Sum[dT, {i, 1, n}], where T0 is the initial temperature of the process and is 500 (kelvin).
4. Δh is an apparent activation energy in J/mol
5. tao,k (from equation 20) is the relaxation time and it is a function of Tf,k-1 (The initial Tf(0)=500) and Tk
6. Tf,k which is what I want to calculate can be considered to be as the "glass transition temperature". Hence, it depends on the material and it is quite different than T,k which is just the process temperature.

So, as you can see, every parameter is known except Tf which is what I want to solve by doing this numerical method given the known parameters.

At time zero, for example

tao = A*Exp[((x*ah)/(R*Tk)) + ((1 - x)*ah/Subscript[T, fkminus1])] = 
   7.6*10^-38*
    Exp[((0.49*342496)/(8.314*499.75)) + ((1 - 0.49)*342496/500)];

and then with the above tao we use it to calculate Tf at time 0 and then you can calculate tao at the next time, which allows you to calculate ‘Tf’ and the same time. The process repeats in an iterative way (which I don’t know how to do in Mathematica) until we are able to plot for instance Tf vs T for a set of T. The set of T in this case goes from 500 to 350 (to be able to plot as in the figure). Notice that because the temperature set is already fixed between 500 and 350, the number of steps n is simply n = (T0 - Tfinal)/dt

ah = 342496;(*J/mol*)R = 8.314;(*J/mol.K*)A = 7.6*10^-38;(*s*)b = 0.67;
x = 0.49;
q = -0.1/60;(*K/s*)T0 = 500;(*K*)Tfinal = 350;(*K*)dt = 300;(*s*)n = 
 IntegerPart[Abs[(T0 - Tfinal)/dt/q]];(*number of steps*)Tf0 = T0;
dT = dt*q;(*K*)T[k_] := T0 + Sum[dT, {i, 1, k}];


tf[0] = T0; Do[
 tau[k] = A*Exp[(x ah)/(R T[k]) + (1 - x) ah/(R tf[k - 1])]; 
 tf[k] = T0 + 
Sum[dT (1 - Exp[-(Sum[dT/(q*tau[j]), {j, i, k}])^b]), {i, 1, k}];, {k, n}]

Tfc = Table[{T[k], tf[k]}, {k, n}]; 
T[0] = T0; cp = 
 Table[{T[k], (tf[k] - tf[k - 1])/(T[k] - T[k - 1])}, {k, n}];

{ListPlot[Tfc, PlotRange -> All, AxesLabel -> {"T", "Tf"}], 
 ListPlot[cp, PlotRange -> All, AxesLabel -> {"T", "cp"}]} 

Tf0 = tf[n];q = 0.1/60;(*K/s*)T0 = 350;(*K*)Tfinal = 500;(*K*)dt = 300;(*s*)n1 = 
 IntegerPart[Abs[(T0 - Tfinal)/dt/q]];(*number of steps*)
dT = dt*q;(*K*)T1[k_] := T0 + Sum[dT, {i, 1, k}]; 

tf1[-1] = Tf0; T1[0] = T0;

tau[0] = A*Exp[(x ah)/(R T0) + (1 - x) ah/(R Tf0)]; Do[
 tau[k] = A*Exp[(x ah)/(R T1[k]) + (1 - x) ah/(R tf1[k - 1])]; 
 tf1[k] = Tf0 + 
   Sum[dT (1 - Exp[-(Sum[dT/(q*tau[j]), {j, 0, i}])^b]), {i, 0, 
     k}];, {k, 0, n1}]

Tfh = Table[{T1[k], tf1[k]}, {k, n1}]; cp1 = 
 Table[{T1[k], (tf1[k] - tf1[k - 1])/(T1[k] - T1[k - 1])}, {k, n1}];

{ListPlot[{Tfc, Tfh}, AxesLabel -> {"T", "Tf"}, PlotLegends -> {"Cooling", "Heating"}, PlotStyle -> {Blue, Red}], 
 ListPlot[{cp, cp1}, AxesLabel -> {"T", "cp"}, 
  PlotLegends -> {"Cooling", "Heating"}, PlotStyle -> {Blue, Red}]}