We solve the five one-dimensional hydrodynamical equations for the evolution of the stellar structure with the mass coordinate as independent variable.
We first convert the five non-linear stellar structure differential equations into difference equations following the Henyey's scheme by defining discrete shells in mass. We then solve the resulting linearized system of equations by the Newton-Raphson relaxation procedure (see e.g. Press et al. 1986). The shell number and thinness in mass is adjusted at each time step in order for the density, temperature, pressure, luminosity, internal energy, nuclear energy production rate and mass fractions profiles to be smooth enough. Typically, 600 shells are required for a model on the main sequence phase, while along the AGB, 1000 to 1300 shells are needed (200 to 500 more shells being added during a thermal pulse).
The time step is constrained in order for the density, temperature or nuclear energy generation to not change by more than 10% in each mass shell between two successive models. In addition, for intermediate phases between central burning phases, the time step cannot exceed the relevant Kelvin-Helmholtz time scale. Typically 200 (700) models are necessary to model the central hydrogen (helium) burning phase. During the AGB phase, 2500 to 4500 models are required to model the time separating two successive thermal pulses, while the thermal pulse itself requires 100 to 350 models. Let us stress that these rather small numbers of mass shells and time steps have been chosen in order to save computer time and present a very large set of nucleosynthesis predictions. Straniero et al. (1996) very recently showed that this is probably in the disadvantage of the occurrence of the third dredge-up along the TP-AGB phase, in particular in the present models.
Finally, the mass 
lost during each time step
is suppressed proportionally to the shell masses of the whole
envelope above the burning regions.
The nucleosynthesis equations are solved using the Wagoner (1969) numerical technique, well suited for our purposes. However, if a diffusion treatment is required to simultaneously compute the time-dependent convective mixing and the nuclear burning, the diffusion Eqs. (4) are also converted into difference equations, as for the stellar structure equations, in order to use the same Newton-Raphson algorithm (a very similar method was already used by Sackmann et al. 1974).