The terms within-individual variance and among-individual variance are not commonly found in the mixed effects model literature. It more commonly arises in the ANOVA literature, and rather than "among", the usual term is "between". Total variance is partitioned into that which is attributable to differences within individuals, for example the natural variation that occurs in the measurement of blood pressure of a person during a day, and that which is attributable to differences between (or among) individuals. Some people have generally different blood pressure than others, but each person's blood pressure also varies throughout the day.
A repeated measures ANOVA can be formulated as a mixed effects model. In the case of repeated measures within individuals, a random intercepts model will estimate a variance at the individual level, which will be the variance of $b$ in your notation above, while the residual variance is at the measurement level, and is the variance of $\varepsilon$ in your model. The former is between (among) and the latter is within.
There are 4 simulations here for each individual: the continuous, time-varying covariate $Z_i(t)$, the time-fixed binary covariate $X_i$, the event time $T_i$ (based on $X_i$, the assumed time-varying coefficient $\phi(t)$, and $Z_i(t)$), and the censoring time $C_i$. The censoring simulation needs to be handled last, after the rest of the simulation. This usually takes some playing around to meet any specific censoring percentage.
Start with simulating $Z_i(t)$ among the individuals, as indicated in the first paragraph of "Theoretical Set up." You need to do that in continuous time for the simulation of event times, but in the final simulated data set for modeling you restrict the recorded values for each individual to the values at the specific observation times $R_t$. If $\max(R_i(t)) \le 20$, don't keep $Z_i(t)$ values in the final simulated data set after the time of the 20th point of the time grid.
As the assumed form of $Z_i(t)$ has period 1 in $t$, I at first assumed that the authors intended t=1 to represent the maximum observation time, $\tau$, the upper limit in the displayed integral. In that case, then the 25 total potential observation times would be seq(from = 0.04, to = 1, length.out=25), and the restriction to at most the first 20 observations would mean that you don't record those discretely sampled values beyond t=0.8 for $Z_i(t)$.
Reconsideration: The proposed baseline hazard function has a median survival of about 5 time units, however, so that initial assumption seems to be wrong. I'm not sure what the authors intended the maximum observation time to be. The basic idea in the prior paragraph still holds, however: you have 25 evenly spaced potential observation times for event occurrence, but you only record values of $Z_i(t)$ for the first 20.
For the integral, once you choose a set of the 3 $\lambda$ values, the integrand has a closed form amenable to numeric integration. As written, it seems that the argument to $\exp$ in the formula for $h(t)$ is fixed, simplifying the subsequent integration of $h(t)$ to get the cumulative hazard $H(t)$.
Reconsideration: The formula for $h(t)$ presented by the authors doesn't make sense for a proportional-hazards model with time-varying covariates and time-varying coefficients. In fact, the formula presented by the authors has a time-fixed covariate for each individual that depends explicitly on the last observation time. I think that they intended to write:
\begin{align*}
h_{i}(t)=h_{0}(t) \exp{\left\{\alpha_{1} X_{i}+ \phi(t) Z_{i}(t) \right\} }
\end{align*}
for the instantaneous hazard, which would then be integrated over time to get the cumulative hazard, up to the last observation time $\tau$.
That handles items 1 and 3, to the limits of what I can figure out. After you have chosen the random $\lambda$ values, $F_i(t)$ for individual $i$ is a continuous closed-form function, used along with the random value of the time-fixed binary covariate $X_i$ to do the integration (in one or the other of the forms discussed above) for each individual's hazard function. But for the $Z_i(t)$ values in the simulated data set, you start by only keeping the values of $Z_i(t)$ up to the 20th discrete observation time.
To simulate each individual's event time, you sample from $U(0,1)$. The authors say to find the corresponding time from the inverse of the individual's cumulative hazard function over time, but I think that they meant to sample from the corresponding survival-time distribution, where $S(t)=\exp(-H(t))$. The integration to get the cumulative hazard only needs to be done out to the last observation time. If the random sampling indicates an event time beyond that, you indicated a censored value at the last observation time.
Only then does it make sense to simulate from $U(0,C_{max})$ for censoring times. There is no way that I know to choose $C_{max}$ simply. Don't forget that some event times will be censored at the last observation time of the study, also. You try some value, find out what fraction of cases would be censored (based on $C_i < T_i$ or $T_i>$ last observation time), and if that doesn't work keep on iterating.
Once you have found a value of $C_{max}$ that gave an appropriate fraction of censored cases, omit from the final data set any data for individual $i$ for which the observation time is greater than the corresponding $C_i$.
Best Answer
The first model says that the hazard of an event at time point $t$ is associated with the underlying level of the longitudinal outcome at $t$. By underlying we mean without the measurement error (biological variation in some contexts). The coefficient $\alpha$ is the corresponding log hazard ratio for one unit increase of the longitudinal outcome at time $t$.
Presumably, in the second model, the term should be $\boldsymbol{\alpha}^\top \mathbf{b}_i$, not $\mathbf{b}_i$. This model says that the hazard of an event depends on the random effects from the model for the longitudinal outcome. The formulation makes more sense when you have only random intercepts or random intercepts and random slopes in the linear mixed model. In this case, you would have $\alpha_1 b_{i0} + \alpha_2 b_{i1}$. The coefficient $\alpha_1$ would be the log hazard ratio for a unit increase in the random intercept $b_{i0}$, and $\alpha_2$ the log hazard ratio for a unit increase in the random slope $b_{i1}$ (given each time that all other covariates and the other random effect are fixed). This formulation loses its nice interpretability when you include nonlinear terms in the design matrix for the random effects $\mathbf{z}_i(t)$ (because then you cannot interpret the $\boldsymbol{\alpha}$'s in isolation).