%\documentstyle[11pt,../par_a4]{../myarticle}
%\documentstyle[referee]{laa}
%\documentstyle{laa}
\documentstyle{l-aa}
%\documentstyle[referee]{l-aa}

\input{psfig.tex}

\def\defbox{\centerline}
%\def\defbox#1{ }

\def\dsdx#1#2{{\partial #1 \over \partial #2}} %d/dx
\def\dsdxx#1#2{{\partial^2 #1 \over \partial #2^2}} %d/dx
\def\ddx#1#2{{d #1 \over d #2}} %d/dx simple
\def\dpar#1#2{{\partial #1/ \partial #2}} 
\def\Etal#1{et al.\ (\cite{#1})}	%avec  parentheses
\def\etal#1{et al.\ \cite{#1}}
\def\dsty{\displaystyle}
\newcommand{\boite}[2]{\centering{\fbox{\rule{#1}{0cm}\rule{0cm}{#2}}}}

\def\mref#1{(\protect\ref{#1})}
\def\fig#1{Fig.\,~\protect\ref{#1}}
\def\Fig#1{Figure\,~\protect\ref{#1}}
\def\tab#1{Table\,~\protect\ref{#1}}
\def\eq#1{Eq.\,~\mref{#1}}
\def\Eq#1{Equation\,~\mref{#1}}
\def\sec#1{Sect.\,~\protect\ref{#1}}
\def\Sec#1{Section\,~\protect\ref{#1}}
\def\app#1{App.\,~\protect\ref{#1}}
\def\App#1{Appendix\,~\protect\ref{#1}}
\def\tot{{\rm tot}}
\def\ext{{\rm ext}}
\def\LoD{{\rm d}}
\def\obs{{\rm obs}}
\def\rmc{{\rm calc}}
\def\bib{\par\noindent\hangindent=3mm\hangafter=1}
\def\Box{\raisebox{.5ex}{\framebox[1ex]{}}}
\def\respt#1{(respt.\ #1)}
\def\eg{e.g.\ }	%exampli grati
\def\etc{etc...\ }
\def\noi{\noindent}
\def\op{op. cit.}	%opus citato
\def\loc{loc. cit.}	%loco citato
\def\ie{i.e.}	%ita est

\begin{document}

\thesaurus{03.13.4 - 06.05.1 -- 06.09.1-- 08.05.3 -- 08.09.3}

\title{CESAM: a code for stellar evolution calculations}
\author{P. Morel}

\institute{Cassini, URA CNRS 1362 and GDR CNRS G 131, Observatoire de la C\^ote 
d'Azur, Nice, France}

\offprints{morel@obs-nice.fr}

%\date{Version: (\today)}

\date{Received date; accepted date
%, {\Large\bf PREPRINT, ONLY SUBMITTED, NOT YET ACCEPTED}
}

\maketitle

\begin{abstract}
The code CESAM is a consistent
set of programs and routines which performs calculations of 1D quasi-static
stellar evolution including diffusion and rotation.
The principal new feature is the solution of the quasi-static equilibrium
by collocation method based on piecewise polynomials approximations projected
on their B-spline basis; that allows stable and 
robust calculations and the exact restitution of the solution not only
at grid points. Another advantages are the monitoring by only one parameter
of the accuracy and its improvement by superconvergence.
An automatic mesh refinement has been designed for adjusting the localizations
of  grid points according to the changes of unknowns, each limit between 
a radiative and a convective zones being shifted to the closest grid point.
For standard models, the evolution of the chemical composition is solved by 
stiffly stable schemes of orders up to four; for non-standard models the 
solution of the diffusion equation employs the Petrov-Galerkin scheme, with
the mixing of chemicals in convective zones performed by strong turbulent
diffusion. A  precise restoration of the atmosphere is allowed for.
CESAM computes evolution of stars from the pre-main sequence 
to the beginning of the $^4$He burning $3\alpha$ cycle.
In this paper a detailed description of the algorithms is presented.

\keywords{Methods: numerical -- sun: evolution -- sun: interior -- stars:
evolution -- stars: interior}
\end{abstract}

\section{Introduction}
The endeavor of a numerical code of stellar evolution is, more or less, 
directed 
towards a given goal; either the purpose is the computation of models evolved 
for long time intervals as for the study of AGB stages (Han \etal{ha}), in which
cases the search 
of efficient and low-cost algorithms is the priority, or a high level of 
accuracy has to be reached, as for solar models which need high numerical
accuracy in addition to 
precise physics (Christensen-Dalsgaard, \cite{cd0}; Reiter \etal{re}).

A first way for improving the numerical precision of a stellar model is the 
refinement of the mesh, a second is the use of a numerical algorithm of high 
order. Until recently the codes used for stellar evolution were, more or less, 
derived from the Henyey's code (Henyey \etal{hl}, \cite{hf}) which uses 
finite differences and so, have a {\em fixed} order of accuracy, \eg first or 
second; with this class of schemes, the only way to 
have various orders of accuracy is to write several algorithms.

For the integration of two-point boundary value problems, high accuracy 
(typically $10^{-10}$ or even better) can be obtained with the multishooting 
technique (Stoer \& Bulirsch, \cite{sb}, Sect. 7.3.5) employed in few 
codes of stellar evolution (Wilson, \cite{wi}; Reiter \etal{re}). They 
are CPU consuming as reported by Embarek (\cite{em}) but they are well designed 
for parallel computers (Reiter et al. \loc). Solar models with such 
impressive 
numerical accuracy have been calculated with a very few number ($\la 50$) of 
fitting points. These methods, however, suffer of a conceptual gap between the 
reachable accuracy in space and in time since the mixing of the 
chemical species in Convective Zones (CZ) with such high-accuracy 
schemes is, in practice, only possible if neither thermonuclear reactions nor 
diffusion take place in CZ; therefore they are well suited for the calculation 
of 
Standard Solar Models (SSM).
Applications of solar models are 
the computation of the eigenfrequencies; that needs a precise knowledge of the 
solution around the nodes 
whose locations change according to the modes, their 
orders and degrees; therefore the solution needs to be restored at locations
well
distributed around the turning points, \ie, between the grid points, that needs
lengthy calculations with  multishooting methods if the {\em full accuracy}
is preserved.

The goal of CESAM\footnote{This acronym means: Code d'Evolution Stellaire 
Adaptatif et Modulaire.} is to perform a standard solar evolution from ZAMS to 
present age with an accuracy of the order of $10^{-4}$, using 
about 500 mass shells and 50 time steps (with a maximum time step of 100 million 
years); according to the physics employed it 
also calculates evolution from PMS to the onset of $^4$He 
burning for various stellar masses. The method of integration does not use the 
familiar finite differences; it is based 
on the so-called splines collocation method derived by de\,Boor (\cite{db}, 
Chapt. XVII). The basic idea is to seek piecewise polynomials as 
solutions of differential problems; among all the bases of the linear space of 
piecewise polynomials, the B-splines basis (de\,Boor, \loc; Schumaker, 
\cite{sk}) has been designed especially for calculation purposes. Its main 
advantages are: (i) the order of the scheme has only one parameter dependence,
(ii) algorithms specially derived for spline calculations allow efficient and 
stable computations, (iii) a convenient method has been designed for 
connecting the pieces of polynomials with a given order of continuity,
(iv) only few changes in the algorithms allow 
to get either the weak solutions (Galerkin's methods) or the
strong solutions (collocation method),
(v) ability of restoring exactly the numerical solution anywhere, (vi) 
ability to handle discontinuous solution, (vii) no need 
of a peculiar algorithm at the center.
However, due to the non-trivial and 
unfamiliar algebra of B-splines, the algorithms are much more intricate than 
with the standard finite differences.
With the use of a functional basis (\ie\ presently, the B-spline basis) it is 
a triviality to structure 
the code in a "numerical space" where the differential problem is formally 
solved and in a "physical space" where the equations are formed whatever the 
numerical method used for their solution is; hence all the physics: Equation Of 
State (EOS), opacities, nuclear reaction rates, external boundary conditions, 
calculation of the convective flux, \etc are calculated in {\em external} 
routines which can be, on need, user supplied or taken from the source itself 
without any change in the "numerical space". Hence CESAM allows calculations of 
stellar models with various physical assumptions, physical data, external 
boundary conditions, numerical methods and numerical accuracy; it also includes
diffusion and mass loss. 



For consistency the same order of accuracy in time and in
space is desirable; with respect to time, the fourth order is a hardly reachable 
upper limit, that is 
due to (i) the stiffness of the problem: the evolutionary 
time scales of the various chemical species differ by more than eighteen orders 
of magnitude, (ii) the convective mixing can hardly be done with 
a high order scheme. For standard models, \ie, without diffusion, the numerical
scheme currently employed is the first order 
Euler's backward formulae (Kippenhahn \etal{ks}), with 
simplifications (Arnett \& Truran, \cite{at}) or improvements (Wagoner, 
\cite{wa}), (Christensen-Dalsgaard, \cite{cd0}) makes use of the second order
trapezoidal rule; CESAM uses a stiffly stable Implicit 
Runge Kutta (IRK) scheme with order up to 4. With diffusion, Runge Kutta scheme
(Alecian, \cite{al}), finite difference schemes (Cox \etal{co}), Crank Nicholson 
(Proffit \& Michaud, \cite{pm}; Charbonel 
\etal{cvz}) are employed; in CESAM 
the diffusion and thermonuclear equations, with mixing of chemicals in CZ, are 
solved as a whole using the implicit
backward first order Euler's scheme in time and, in space,
the Petrov-Galerkin's method with order up to four.

Some recent papers are based on calculations made with CESAM (Morel \etal{mp1}; 
Morel \etal{mp2}, \cite{mp4}; Berthomieu \etal{be}; Goupil \etal{go}; Audard \& 
Provost, \cite{ap}; Lebreton \etal{le}; 
Chmielewski \etal{ch}; Morel \& Schatzman, \cite{mp6}). Preliminary and short 
descriptions of the numerical methods have already been given (Morel, 
\cite{mp0}; Morel, \cite{mp5}, \cite{mp3}) or referenced (Gabriel, \cite{gab}); 
the algorithms have been 
adapted for the computation of evolution of giant planets (Guillot \& Morel, 
\cite{gm}). This 
paper is devoted to a detailed description of the numeric of CESAM.


In \Sec{sec:pd} the basic physics of stellar evolution is briefly recalled for 
references and notations; special emphasis is given to the mixing of CZ, on 
initial 
conditions and on external boundary conditions 
for which a detailed restoration of an atmosphere is allowed for. \Sec{sec:ivp} 
is devoted to the numerical techniques; 
a first part is concerned with the solution of the two points boundary problem
with emphasis on the choice of the integration variables and on the automatic 
location of grid points; in a second part, the schemes used for solving 
the 
stiff differential initial value problem of the evolution of the chemical
composition with, and without, diffusion are 
described. In \Sec{sec:gf} some indications concerning the management of the code 
are given. Results of investigations on internal accuracy are presented
in \Sec{sec:test} and the discussion in \Sec{sec:dis}.
Three appendices are devoted to numerical technical discussions.


\subsection{Contents of the paper}
\begin{itemize}
 \item \sec{sec:pd}, Physical description
  \begin{itemize}
   \item \sec{sec:es}, The equations for stellar evolution
   \item \sec{sec:dc}, Diffusion of chemicals
   \item \sec{sec:conv}, Convection
   \begin{itemize} 
    \item \sec{sec:mix}, Mixing
    \item \sec{sec:ovsht}, Overshooting, undershooting 
   \end{itemize}
   \item \sec{sec:bc}, Boundary conditions
   \begin{itemize}
    \item \sec{sec:ml}, Mass loss and mass change
    \item \sec{sec:bcc}, Boundary conditions at center
    \item \sec{sec:ebc}, External boundary conditions
   \end{itemize}
   \item \sec{sec:ic}, Initial conditions  
   \begin{itemize}
    \item \sec{sec:pms}, PMS initial models
    \item \sec{sec:zams}, ZAMS initial models
   \end{itemize}   
  \end{itemize}
  \item \sec{sec:ivp}, Solving the two points boundary initial value problem
  \begin{itemize}
   \item \sec{sec:oni}, Overview on numerical integration
   \item \sec{sec:qse}, Numerical solution of the quasi-static equilibrium
   \begin{itemize}
    \item \sec{sec:lv}, The lagrangian variables employed
    \item \sec{sec:agp}, Automatic location of grid points
    \item \sec{sec:bccn}, Boundary conditions at center
    \item \sec{sec:ebcn}, External boundary conditions
    \item \sec{sec:tds}, Gravothermal energy
    \item \sec{sec:eul}, An eulerian set of variables
   \end{itemize}
   \item \sec{sec:evol}, Numerical solution of the evolution of chemicals
   \begin{itemize}
    \item \sec{sec:chim}, Interpolation of chemicals
    \item \sec{sec:mixc}, Standard evolution
    \item \sec{sec:mixd}, Evolution with diffusion      
   \end{itemize}
  \end{itemize}
  \item \sec{sec:gf}, General feature for the implementation of the code
  \begin{itemize}
   \item \sec{sec:upd}, Initialing and updating the solution
   \item \sec{sec:tsc}, Time step controls
   \item \sec{sec:tv}, Sets of numerical parameters
   \item \sec{sec:cal}, Calibration of solar models
  \end{itemize}
  \item \sec{sec:test}, Estimates of the numerical accuracy
  \begin{itemize}
   \item \sec{sec:ctes}, Physics used for the numerical tests 
   \item \sec{sec:scv}, Superconvergence
   \item \sec{sec:iacc}, Estimate of the internal accuracy
   \item \sec{sec:dcs}, Discontinuities
   \item \sec{sec:cev}, Comparison of evolutionary tracks
  \end{itemize}  
  \item \sec{sec:dis}, Discussion
  \item \app{sec:summ},  Summary of abbreviations and notations
  \item \app{sec:ap1}, Sets of variables without central singularities
  \item \app{sec:set}, Setting a grid point at a given location
  \item \app{sec:pbc}, Precise external boundary conditions
  \item \app{sec:sirk}, Stable integration with IRK formulas 
% \item \app{sec:imx}, Instantaneous mixing with IRK formulas
  \item \app{sec:ided}, Integration of the diffusion equation
  \item \app{sec:ap2}, Calculations with B-splines
  \begin{itemize}  
   \item \app{sec:nod}, The nodal vector: the key of the computations with 
B-splines
   \item \app{sec:int}, Interpolation
   \item \app{sec:bsc}, Computation of the B-spline coefficients
   \item \app{sec:cpw}, Piecewise polynomial interpolation is conservative 
   \item \app{sec:ide}, Integration of differential equations
   \begin{itemize} 
    \item \app{sec:dbb}, Collocation with the de\,Boor's basis
    \item \app{sec:ga}, Basis for Petrov-Galerkin's method
    \item \app{sec:sol}, Solving the linearized systems
   \end{itemize}     
  \end{itemize} 
    
\end{itemize}

\section{Physical description}\label{sec:pd}
\subsection{The equations for stellar evolution}\label{sec:es}
Although several physical laws and data are 
more easily derived from the set of variables of state $\{ \rho,T,{\cal X} \}$ 
(\ie
density, temperature, vector of the abundances of chemical species), with the 
numerical techniques employed, the set $(P,T,{\cal X})$ is more 
suitable\footnote{Recall that, due to the mixing of chemicals a discontinuity 
may occur in the density profile at the limit of convective cores.}; 
$P$ is the pressure.

Let, as independent lagrangian variable, $M$ be the mass inside the sphere of 
radius $R$; at first order for the rotation, in the assumption of spherical 
symmetry, without magnetic field, the system of 1D 
rotating stellar structure equations:
\begin{equation}
\label{eq:eq1}\displaystyle
 \left\{
 \begin{array}{l}\vspace{0.1 in}
  \dsdx P M = -\frac {GM}{4\pi R^4} + \frac{\Omega^2}{6\pi R} \\ \vspace{.1 in}
  \dsdx T M =  \dsdx P M \frac T P \nabla \\ \vspace{.1 in}
  \dsdx R M =  \frac 1{4\pi R^2\rho} \\ \vspace{.1in}
  \dsdx L M =  \epsilon -\dsdx U t +\frac P{\rho^2}\dsdx\rho t\\ \vspace{.1in}
  \dsdx {X_i}t = -\dsdx {F_i}M + \Psi_i,\ 1\leq i\leq n_{\rm X}
  \end{array}
 \right.
\end{equation}
\noi are recalled for notation and reference. Standard symbols 
are employed; $t$: time, $G$: gravitational constant, $\Omega(M,t)$: angular 
velocity, 
$\nabla\equiv\dpar{\ln T}{\ln P}$: temperature gradient, $L(M,t)$: luminosity, 
$\epsilon$:
rate of nuclear energy release, $U$: internal energy per unit of mass, 
$X_i(M,t)$: fraction of 
unit mass which consists of nuclei of type $i$, ${\cal X}\equiv 
(X_1,\ldots,X_{n_{\rm X}})^T$, $F_i(M,t)$: flow of $X_i$ due to 
diffusion, $\Psi_i$: creation rate of $X_i$ due to thermonuclear reactions, 
$n_{\rm X}$: number of chemical species.

\subsection{Diffusion of chemicals}\label{sec:dc}
For the chemical species $i=1,\ldots,n_X$, in its most general form, the flow 
$F_i(M,t)$ is written:
% (Proffit \& Michaux, \cite{pf})
\begin{equation}\label{eq:fld}
 F_i=-4\pi R^2\rho\left( 4\pi R^2\rho D_i\bullet\nabla\!_{\rm M}{\cal 
X}-v_iX_i\right),
\end{equation}
\noi the components $d_{i,j},\ i\not= j,$ of the vector: 
\[
D_i(M,t)\equiv (d_{i,1},\ldots,d_{i,n_X})^T,
\] 
\noi are the coefficients, $d^\star_{i,j}$, of microscopic atomic 
diffusion for the chemical species \#$i$ with 
respect to the species \#$j$; the $(i,i)$-th component: 
\begin{equation}\label{eq:d}
d_{i,i}\equiv d^\star_{i,i}+d_{\rm T}
\end{equation}
\noi also includes the turbulent diffusion $d_{\rm T}$; $v_i(M,t)$ 
is the velocity of diffusion;
the "$\bullet$" denotes the inner product:
\[
D_i\bullet\nabla\!_{\rm M}{\cal X}\equiv \sum_j d_{i,j}\dsdx{X_j}M.
\]
\noi  Recall that, owing to diffusion, the chemical 
composition is a continuous function.

\paragraph{Extended vector of chemicals.} 
With the diffusion of $Z$, the content of heavy elements per unit of mass, 
and/or of ${\cal 
M}_\Omega\equiv\Omega R^2$, the angular momentum per unit of mass, it is 
convenient to generalize the vector of chemicals by adding one/two 
components, since, with {\it ad hoc} coefficients, it fulfills a diffusion 
equation similar to \eq{eq:d}. Hereafter, otherwise stated, for sake of 
simplicity, the 
previous notations, \ie, $\cal X$ and $n_X$, will be also used, respectively, for 
the extended vector of chemicals and for its length:
\[
{\cal X}={\cal X}(M,t)\equiv (X_1,\ldots,X_{n_{\rm X}},Z,{\cal M}_\Omega)^T.
\]
\subsection{Convection}\label{sec:conv}
At each level, according to a choice made {\it a priori},
either the Schwarzchild's or the Ledoux's criterion allows 
to decide if the energy transport is made by radiation or by convection.

\subsubsection{Mixing}\label{sec:mix}
In convective zones, contingently extended by overshooting (hereafter named 
"Mixed Zones", MZ), the characteristic time scale of the convection is small 
with respect to the thermonuclear time scales, then the MZ are {\em quasi 
instantaneously} mixed and their chemical composition is homogeneous. Owing to 
thermonuclear reactions, at each Limit between a Mixed zone and a Radiative zone 
(LMR), the chemical composition may have discontinuity. It is emphasized that, 
if it is not regenerated, a long-lasting discontinuity on 
chemicals is an unphysical situation.
With CESAM, the mixing is made either by standard integration, 
\ie, instantaneous mixing, or by diffusion.

\paragraph{Instantaneous mixing.}\label{sec:mixisd}
For a given MZ, let $F_i^+$ and $F_i^-$ be the flows of $X_i$ at each LMR, the 
last equation of the system 
\eq{eq:eq1} 
becomes (Clayton, Sect. 6.3, \cite{cy}):
\begin{equation}\label{eq:dx}
 \ddx {\bar{X_i}}t=\frac{F_i^++\int_{\rm MZ}\Psi_i(P,T,\bar{\cal X},t) {\rm d}M
  +F_i^-}{\int_{\rm MZ} {\rm d}M}
\end{equation}
\noi where each component $\bar X_i$ of $\bar{\cal X}$, the vector of the mean 
abundances, is expressed as:
\begin{displaymath}
 \bar X_i=\int_{MZ}X_i {\rm d}M\ / \ \int_{MZ} {\rm d}M.
\end{displaymath}

\noindent The reaction rates, $\Psi_i(P,T,{\cal X},t)$, are non-linear functions 
of the abundances, therefore, at any point $M$ of a MZ, the rate $\Psi_i$ 
has to be calculated with $\bar{\cal X}$ the {\em mean values, at time $t$,} of 
the 
abundances over the MZ; hence the calculation of the chemical composition is 
an integro-differential problem.

\paragraph{Mixing by diffusion.}\label{sec:mixiad}
As suggested by Eggleton (\cite{eg}, \cite{eg1}), the convective mixing can be 
modelled by {\em turbulent} diffusion; that can be made by adding a "mixing" 
diffusion coefficient $d_{\rm M}(M,t)$ to the turbulent diffusion $d_{\rm T}$; 
\eq{eq:d}  becomes:
\begin{equation}\label{eq:d1}
d_{i,i}\equiv d^\star_{i,i}+d_{\rm T}+d_{\rm M}.
\end{equation}

\paragraph{Mixing of angular momentum.}
In MZ solid rotation is assumed, there, the angular velocity is written:
\begin{displaymath}
\Omega(m,t)=\frac{\int_{{\rm MZ}(t)} r^2(m,t-\Delta t)\Omega(m,t-\Delta t)
{\rm d}m}{\int_{{\rm MZ}(t)} r^2(m,t){\rm d}m}.
\end{displaymath}
\noi with the integrals taken over all the MZ at time $t$; $\Delta t$ is the 
time step.
Likewise, when the diffusion of the angular 
momentum is taken into account, ${\cal M}_\Omega$ is homogenized, in each MZ, 
by strong turbulent diffusion. 

\subsubsection{Overshooting, undershooting}\label{sec:ovsht}
At LMR the velocity of 
convective eddies is not zero, therefore the convective motions can penetrate 
into 
the radiative zone. Owing to the small characteristic time scale of the 
convection, the matter is homogenized within the overshoot/undershoot zone;
if, there, the temperature gradient keeps its convective adiabatic value with
a discontinuity at the limit of the extension,
one speaks of "convective penetration" (Zahn, \cite{jpz}).

\subsection{Boundary conditions}\label{sec:bc}
With respect to the independent variables $(M,t)$, boundary and initial
conditions for the dependent variables $P$, $T$, $R$ and $L$ must be specified.


\subsubsection{Mass loss and mass change}\label{sec:ml}
Let $\dot M\leq 0$ be the rate of external mass loss due to stellar wind and
$\dot{\cal M}\leq 0$ be the change of mass due to mass defect (for the Sun 
$\dot{\cal
M}\simeq L/{\cal C}^2\sim 10^{-14}M_\odot y^{-1}\sim \dot M$):
\begin{displaymath}
\dot{\cal M}(t)=\int_0^{M_\star(t)}\frac{\epsilon+\epsilon_\nu}{{\cal C}^2}{\rm 
d}m
\end{displaymath}
\noi $\epsilon_\nu$ is the energy of neutrinos, $\cal C$ the velocity of light
and $M_\star(t)$ the total mass of the star at time $t$:
\[
M_\star(t)=M_{\rm ini}+\int_0^t\left(\dot M(x)+\dot{\cal M}(x)\right) dx
\]
\noi with $M_{\rm ini}$ as the initial mass of the star.
It is assumed that the mass lost by stellar wind is simply detached from the 
star by some process which is not described.

\paragraph{Loss of angular momentum.} Likewise, per unit of time, the amount of 
angular momentum lost, caused only by mass loss,
amounts to:
\[
\dot{\cal M}(t)_\Omega =\dot M(t) \Omega(M_\star,t)R_{\rm A}^2(t),
\]
with $R_{\rm A}$ as the Alf\`en radius, where matter decouples from the magnetic 
field (Schatzman, \cite{s1}). 


\subsubsection{Boundary conditions at center}\label{sec:bcc}
Due to the hypothesis of spherical symmetry, the radius, luminosity
and the fluxes of chemical and of angular momentum vanish at the 
center, \ie, at $M\equiv0$:
\begin{equation}\label{eq:bc}
 R(0,t)=0,\ L(0,t)=0,
\end{equation}
\noi and:
\begin{equation}\label{eq:bce}
 F_i(0,t)=0,\ i=1,\ldots,n_{\rm X},\  F_\Omega(0,t)=0.
\end{equation}

\subsubsection{External boundary conditions}\label{sec:ebc}
The problem of the external boundary conditions is more intricate. Inside a 
star, in the optically thick medium, the radiative flux can be derived from the 
diffusion equation which is the limit of the transfer equation for large optical 
depth (Mihalas, \cite{ml}, Sect. 2.5). Since it is beyond our possibilities to 
include a full atmosphere calculation into a whole stellar model, one separates 
the internal structure calculations where the diffusion ap\-pro\-xi\-ma\-tion
is valid, from an outer part where the atmosphere calculation is made under 
simplified assumptions; in that sense the atmosphere is only "reconstituted". 
The limit between
the two parts, where the bottom of the atmosphere coincides with the top of the
envelope, must be located at a Rosseland optical depth $\tau_{\rm b}$ large 
enough for the 
diffusion approximation to be valid.
Precise radiative transfer calculations (Morel \etal{mp2}; \cite{mp4}) 
have shown that $\tau_{\rm b}\geq 20$
if one needs to be sure that the diffusion approximation strictly satisfied
at all wavelengths; therefore, owing to convection, the use of purely radiative 
${\cal
T}(\tau)$ laws, such as Eddington's law (Maeder \& Meynet, \cite {mm}), are
not valid for stars of  mass $M_\star\lse 1.8M_\odot$ which have always a 
significant external convective zone. 

At time $t$, the bottom of the atmosphere is connected at the 
outermost boundary of the envelope by three functions of luminosity, radius and
time: pressure $P_{\rm b}(R,L,t)$, temperature $T_{\rm b}(R,L,t)$ and mass 
$M_{\rm b}(R,L,t)$. These three functions are provided by the reconstitution of 
the atmosphere. At any point in the atmosphere,
the chemical composition is assumed to be the chemical composition of the outer
layer of the envelope.
It is emphasized that the differential problem \eq{eq:eq1} has, then, an 
external
{\em open limit}, since $M_{\rm b}$ changes with time though it remains very 
close to the total mass of the star, $M_\star(t)$, typically $M_\star(t)-M_{\rm
b}(R,L,t)\sim 10^{-10} M_\star(t)$; a very 
elegant solution to the numerical problem of the open limit is obtained via the 
automatic location of grid points as seen \sec{sec:aloc} . 


\paragraph{Simplified external boundary conditions.}
For some applications only an estimate of the radiative energy transfer in the 
outer layers is sufficient and good results are obtained with the single shell 
approximation as described by Christensen-Dalsgaard (\cite{cd1}). Working 
with the equation of hydrostatic equilibrium, the definition of the Rosseland 
optical depth, ${\rm d}\tau=-\kappa\rho{\rm d}R$ with $\kappa$ as the Rosseland
mean opacity, setting the definition of the effective temperature $T_{\rm 
eff}$ as $T_{\rm b}$, the simplified external boundary
conditions are written:
\begin{equation}\label{eq:sext}\displaystyle
\begin{array}{l}\vspace{0.1 in}
 \ddx P\tau=\frac{GM_{\rm b}}{\kappa_{\rm b} R_{\rm b}^2}-
 \frac{2\Omega^2 R_{\rm b}}{3\kappa}
\sim\frac{P_{\rm ext}-P_{\rm b}}{\tau_{\rm ext}-\tau_{\rm b}}
=\frac{P_{\rm b}}{2/3}, \\
\ T_{\rm b}=T_{\rm eff}=4\pi\sigma R_{\rm b}^2L_{\rm b},
\end{array}
\end{equation}

\noi here it is assumed that the pressure cancels, $P_{\rm ext}=0$, at the 
optical depth, $\tau_{\rm ext}=0$, of the outer limit of the atmosphere,  
likewise the optical depth of the bottom of the atmosphere is assumed to be 
$\tau_{\rm b}=2/3$;
$\sigma$ is the Stefan-Boltzman constant. The mass $M_\star$, the radius 
$R_\star$ and the luminosity $L_\star$ of the star respectively then is 
written: 
\begin{equation}\label{eq:star}
 M_\star=M_{\rm b},\ R_\star=R_{\rm b},\ L_\star=L_{\rm b}.
\end{equation}
The unknown, $P_{\rm b}(t)$, $T_{\rm b}(t)$ and $M_{\rm b}(t)$ are solutions 
of the non-linear set of \eq{eq:sext}, owing to the conditions \eq{eq:star}.

\paragraph{Precise external boundary conditions.}
For other applications, as for helioseismology purposes, since the  eigenmodes 
are reflected in the atmosphere, for each level, one needs accurate values of 
the physical quantities of necessity for the calculation of the frequencies, \ie, 
a detailed "reconstitution" of the atmosphere is required;
a standard technique is the use of a 
so-called ${\cal T}(\tau,T_{\rm eff},g)$ law (Henyey \etal{hl}; 
Kippenhahn \etal{ks}); here, as far as the atmosphere has a small extent, the 
gravity $g$ is constant within and equal to the total gravity 
of the star, $g_\star\equiv GM_\star/R_\star^2$.

A difficulty with the detailed reconstitution of the atmosphere is connected to
the location of the stellar radius $R_\star$ which encloses the mass of the 
star
$M_\star$; using the standard definition (Schatzman \& Praderie, \cite{sp}), at 
time $t$, the radius of the star is defined at the level of optical depth 
$\tau_\star$ where the temperature 
is equal to the effective temperature:
\begin{equation}\label{eq:star1}
{\cal T}(\tau_\star,T_{\rm eff},g_\star)=T_{\rm 
eff}=\sqrt[4]{\frac{L_\star}{4\pi\sigma R^2_\star}}.
\end{equation}
According to model atmospheres calculated with the Kurucz's ATLAS9 code 
(Kurucz, \cite{kz}), for a 1$M_\odot$  model, from Pre-Main Sequence (PMS) 
initial conditions to present age, $\tau_\star$ varies within the interval
$\sim[2/5~,~2/3]$, since $T_{\rm eff}$ and $g_\star$ vary. 
In some cases, \eg along the evolution of the Sun from the ZAMS to present time, 
the changes of $\tau_\star$
remain small and can be neglected, but it is not the case if PMS 
or post main-sequence are involved -- see Morel \Etal{mp4}, for a 
more complete discussion.

With a time dependent stellar radius, the set of differential equations to be 
solved for the reconstitution of the atmosphere is written:
\begin{equation}
\label{eq:sysa}
\left\{
\begin{array}{l}\vspace{0.05 in}
\ddx P\tau  = \frac{GM}{\kappa R^2} -\frac{2\Omega^2 R}{3\kappa} \\ \vspace{0.05 
in}
\ddx R\tau  = -1\,/\,{\kappa\rho} \\ \vspace{0.05 in}
\ddx M\tau  = -{4\pi R^2}\,/\,\kappa \\ \vspace{0.05 in}
\ddx {R_\star}\tau =  0 \\ \vspace{0.05 in}
\ddx {\tau_\star}\tau  =  0 \\
\end{array}
\right.\left\{
\begin{array}{l}
{\rm boundary\ conditions} \\ \vspace{0.05 in}
\rho(\tau_{\rm min})  = \rho_{\rm ext} \\ \vspace{0.05 in}
M(\tau_\star)  =  M_\star \\ \vspace{0.05 in}
R(\tau_\star)  =  R_\star \\ \vspace{0.05 in}
\vspace{0.05 in}
T_{\rm eff}  =  {\cal T}(\tau_\star,T_{\rm eff},g) \\ \vspace{0.05 in}
R(\tau_{\rm b})  =  R_{\rm b} \\
\end{array}
\right.
\end{equation}
Following the standard method (Stoer \& Bulirsch, \cite{bs}, \S 7.3.0), the two 
last differential 
equations have been added for solving the differential problem with the open 
inner boundary at 
$\tau=\tau_\star$; they allow to link $R_\star$ and $\tau_\star$ to the other 
unknowns. 
At the optically thin limit of the atmosphere, defined by 
$\tau=\tau_{\rm min}$, the density has the value $\rho_{\rm ext}$ 
given by the model atmosphere; here one assumes
$\tau_{\rm min}\equiv 10^{-4}$, that roughly corresponds to the optical depth
of the temperature minimum of the Sun photosphere; hence the model does not 
include the non LTE outerparts of the atmosphere where the density and 
temperature values take place outside of the ranges of the standard opacity and 
EOS tables.
Boundary conditions are therefore located at three different
levels: $\tau=\tau_{\rm b}$, $\tau=\tau_\star(t)$ and $\tau=\tau_{\rm min}$; 
indeed, $\rho_{\rm ext}$ is a function of time since $\rho_{\rm 
ext} (g_\star(t))$.
At time $t$, the solution of \eq{eq:sysa} gives $P_{\rm b}$, $T_{\rm b}$ and
$M_{\rm b}$ for given $R_{\rm b}$ and $L_{\rm b}$.

Obviously, the $T(\tau)$ laws and the internal structure need to be computed 
with consistent physics (EOS, opacity, convection theory, mixing-length 
parameter, chemical composition \etc).


\paragraph{External boundary conditions for diffusion of chemicals.}
As far as the amount of material in the atmosphere remains negligible the 
changes of abundances caused by diffusion can be ignored and a homogeneous 
chemical 
composition can be assumed there. As already stated, the mass 
loss by mass lost and also the angular momentum lost, are here assumed
simply detached by processes which are not described, 
at the level $M=M_{\rm b}(t)$, due to the mass conservation, there is no input 
nor output of chemicals and angular momentum therefore: 
\begin{equation}\label{eq:bcn}
F_i(M_{\rm b}(t),t)=0,\ i=1,\ldots,n_{\rm X}.
\end{equation}

 
\subsection{Initial conditions}\label{sec:ic}
At time $t\equiv 0$, one has to specify the values of pressure, $P(M,0)$, 
temperature $T(M,0)$, radius $R(M,0)$, luminosity $L(M,0)$ and vector of 
chemical ${\cal X}(M,0)$, at any location within the
star interior, $M\in[0,M_{\rm b}(t=0)]$; they 
are derived from an 
homogeneous stellar model of, either, PMS, or, Zero Age 
Main-Sequence (ZAMS).

\subsubsection{PMS initial models}\label{sec:pms}
Along the PMS evolution the radiative losses are only sustained by the release 
of gravothermal energy; 
as the chemicals have no time dependence (Reiter \etal{re1}), the gravothermal 
energy is 
proportional to the partial derivative of the entropy with respect to time:
\[
\epsilon_{\rm G}\sim-T\dsdx St.
\]
\noi As a consequence, for a zero age PMS model, though in 
quasi-hydrostatic equilibrium, the gravothermal energy cannot 
be taken as 0, as it does for ZAMS models. At the onset of the quasi-static
PMS phase, \ie, at the beginning of the Hayashi line (Kippenhahn \& 
Weigert, \cite{kw}), the stars are fully  convective except, perhaps, in a small 
outer region (Iben, \cite{ib}); the entropy is then constant within the star,
therefore the equation of energy can be written:
\begin{displaymath}\label{eq:ct}
 \dsdx LM=\epsilon_{\rm G}\sim-T\dsdx St=cT
\end{displaymath}
\noi here $c(t)$ is the so-called contraction factor; it characterizes the fully 
convective model in quasi-static equilibrium; for stars with the empirical value 
$c=0.02L_\odot M_\odot^{-1}$K$^{-1}$ the starting point of PMS evolution is 
located in the upper right corner of the HR diagram; it can be adjusted to its 
requested place just by changing $c$.

For a given star, let an 
initial PMS model \# 1 be characterized by a contraction factor $c_1$ and a 
model \# 2 with $c_2$ ($c_2\not= c_1$); between models 1 and 2, the change of 
gravitational energy is equal to the radiative losses during the time step, 
$\Delta 
t$, which, approximatively, amounts to:
\begin{equation}\label{eq:dt}
 \Delta t\sim2GM^2\frac{|R_1-R_2|}{(L_1+L_2)R_1R_2}.
\end{equation}
\noi Therefore \eq{eq:dt} gives an approximate value for the evolutionary time 
scale $\Delta t$ of the model characterized by the contraction factor $c\simeq 
(c_1+c_2)/2$; typically one uses: $c_2\sim 1.1c_1$.

\paragraph{Transitions between pre main-sequence, main-sequence and post 
main-sequence.} 
It is assumed that a PMS model becomes a main-sequence model as soon as 99\% of 
the energy generated has a nuclear origin; at that time,
the model is chemically inhomogeneous, in opposition ZAMS models are
homogeneous. Likewise, as soon as, at center, the abundance of
$^1$H$\leq 10^{-5}$ a main-sequence model becomes a POST main-sequence model.

\subsubsection{ZAMS initial models}\label{sec:zams}
ZAMS models are homogeneous and  in hydrostatic equilibrium, \ie, they fulfill 
the equations \eq{eq:eq1} with the gravothermal 
term equal to zero and homogeneous chemical composition. A ZAMS model is not a 
physical situation because, during the interval of time elapsed from the
instant of the ignition of the thermonuclear reactions which 
occurs when the central temperature reaches $T_c\sim 10^6$K, to the instant
when physical conditions similar to the ZAMS are established -- \ie, when
$T_c\sim 15\,10^6$K, for $1M_\odot$ -- the chemicals 
was not processed by amounts corresponding to the values 
of density and temperature (it has been found (Rossignol, \cite{ro}) that, 
everything else fixed, for ages larger than $\ga 150$My, the internal 
structures of a solar model evolved from ZAMS and of a solar model evolved from 
PMS but older by $\sim30$My, are very similar).




\section{Solving the two-point boundary initial value problem}\label{sec:ivp}
The equations of stellar structure, \eq{eq:eq1}, and their 
numerical counterpart, \eq{eq:syq}, are non-linear; in CESAM they are solved
using the so-called {\em damped Newton-Raphson} method, which 
differs from the standard relaxation method (Clayton, \cite{cy}, chap. 6.4) by 
the fact that the change in the unknowns 
are reduced; that allows to follow the "good descent direction" given by the 
jacobian 
of the linearized equations (Conte \& de\,Boor, \cite{cdb}, \S 5.2) of the 
non-linear problem solved as an optimization problem; doing 
that the time step remains of reasonable size, even with shell sources.

\subsection{Overview on numerical integration}\label{sec:oni}
With instantaneous mixing, the set of partial differential equation 
\eq{eq:eq1} is an initial boundary value, integro-differential
problem\footnote{Formally, with diffusion and turbulent mixing, the 
integro-differential character disappears, the problem is only 
differential.}; for the numerical solution, following Henyey \Etal{hl}, it is 
split into two differential problems:

\begin{enumerate}
 \item  The "evolution of chemicals" in this first step, the chemical 
composition is 
updated for the time step $\Delta t$, taking into account the mixing 
in MZ; in this first step, only $\cal X$ changes, the other variables, \ie, 
$P,T,R,L$, remain fixed. 

 \item The "quasi-static equilibrium", in this second step, the dependent 
variables, \ie, pressure, temperature, radius and luminosity, are adjusted to 
fulfill the quasi-static equilibrium with the value of the chemical 
composition, $\cal X$, previously updated.
\end{enumerate}
For each time step $\Delta t$, starting from an approximate solution, 
the two problems are solved sequentially until global convergence; the whole 
process is fully implicit; {\em all} the variables $P$, $T$, $R$, $L$, $\cal 
X$ are centered at time $t$ (Bressan \Etal{br} center the chemical composition
at half time step, that save computing effort but can generate instabilities
in case of stiffness).
%Note also that, as pointed out by Reiter \Etal{re1}, the convergence 
%does not mean that the accuracy of the closure is reached for any variable
%whatever $M\in[0,M_{\rm b}]$ is.

\subsection{Numerical solution of the "quasi-static equilibrium"}\label{sec:qse}
The numerical solution of the quasi-static equilibrium is made using the 
B-spline collocation method described in \App{sec:ide}; the basic idea is to 
seek 
the unknown functions as piecewise polynomials which are, by turn, projected
on their B-spline basis of order $h_q$ (\ie, of degree $h_q-1$), therefore the
(non-linear) discretized 
differential equations only involve the projections and are solved for them;
finally, their knowledge allows to rebuilt the solution, \ie, the piecewise
polynomials, for any location.

\subsubsection{The lagrangian variables employed}\label{sec:lv}
With the mass as the independent variable, the three first equations of the
differential system \eq{eq:eq1} have a singularity at center, $M\equiv 0$; 
to overcome the difficulty, Taylor series have been employed (Kippenhahn 
\etal{ks}; Paczynski, \cite{pa}). Here an alternate approach is applied:
first, variables which avoid the central singularities are employed and, 
second, the equations are not written right at center, as quoted 
\App{sec:ide}; 
therefore, for the first shell, there is no need of a peculiar algorithm.
Owing to the exponents 2/3 and 2 which are optimal (see \App{sec:ap1}),
among all sets of variables which avoid the central singularities, the 
Eggleton (\cite{eg}) variables\footnote{If negative values of luminosity are 
expected, the 
eulerian set of variables given \sec{sec:eul}, ought to be employed.}:
\begin{equation}\label{eq:eq3}
 \left(\frac M{M_\odot }\right)^{\frac 23},
 \ \left(\frac R{R_\odot }\right)^2,
 \ \left(\frac L{L_\odot }\right)^{\frac 23},
\end{equation}
$L\geq 0$, lead to the most precise numerical solution. From the surface to the
center, the pressure and the temperature typically change,
respectively, by eighteen and five orders of magnitude, therefore logarithms
are used for those variables. {\em In fine,} the following set of lagrangian
variables, free from central singularities, is employed: 

\begin{equation}\label{eq:var1}
\begin{array}{l} 
\xi=\ln P,\ \eta =\ln T,\\
\ \mu=\left(\frac M{M_\odot }\right)^{\frac 23},
\zeta =\left(\frac R{R_\odot }\right)^2,
\ \lambda =\left(\frac L{L_\odot }\right)^{\frac 23}.
\end{array}
\end{equation}

\subsubsection{Automatic location of grid points}\label{sec:agp}
\label{sec:aloc}
Due to large gradients the density of mesh points cannot be kept constant along
the model; moreover, in the course of an evolution, the areas affected by large 
gradients 
change, therefore the grid cannot kept fixed.
The most usual method for grid refinement consists by adding or subtracting 
points according to criterions (see, for example, Henyey \etal{hl}, 
Kippenhahn \etal{ks}). Though this method is very simple in its principle, 
in practice, it is difficult to be properly implemented; an automatic mesh 
refinement is more convenient.
Here, at a given time $t$, the number $n(t)$ of mesh points is given and their 
locations are fixed by fulfilling the condition that, from one grid point to the 
next, the jump of a {\em 
strictly
monotonous} "distribution function", $Q(\mu,t)$, is equal to a "distribution 
constant" $C(t)$ (Eggleton, \cite{eg}; Press \etal{pf}, sect. 16.5); at each 
time $t$ the distribution of the grid points, $\mu_i,i=1,\ldots ,n$, therefore 
satisfies:

\begin{equation}\label{eq:eqq}
Q(\mu_{i+1},t)-Q(\mu_i,t)\equiv C(t),\quad i=1,\ldots ,n-1.
\end{equation}
The choice of $Q(\mu,t)$ is based on an {\it a priori} knowledge of the
behavior of the solution. For each value $t$ of the time, one defines an 
"index" function  $q(\mu,t)$ mapping 
$[1,\mu_{\rm b}]$ on $[1,n]$; the shell with the index 1 \respt{$n$}
corresponds to the center \respt{to the top of the envelope}.
Therefore the integration is made on an {\em equidistant grid}.
In terms of the derivative of $Q$ with respect to $q$:

\begin{equation}\label{eq:psi}
\psi (t) \equiv \left( \dsdx Q q \right)_t,
\end{equation}
\eq{eq:eqq} becomes:

\begin{equation}\label{eq:dpsi}
\left( \dsdx \psi q \right)_t=0.
\end{equation}
Note that $Q$ is a linear function of $q$ as soon as the
Eqs.~(\ref{eq:psi}) and (\ref{eq:dpsi}) are fulfilled.
The change of variables $\mu\rightarrow q(\mu,t)$ gives:
\begin{displaymath}
 \psi(t)=\theta\left(\dsdx\mu q\right)_t,
\end{displaymath}
where: 
\begin{equation}\label{eq:theta}
 \theta(\mu,t)\equiv\left(\dsdx Q\mu\right)_t,
\end{equation}
can be derived from the analytic form of $Q(\mu,t)$.
Thus, there are two more unknowns:
$\psi(t)$ and $\mu(q,t)$; they fulfill a system
of differential equations of first order with boundary conditions:

\begin{equation}\label{eq:muq}
\left(\dsdx \mu q \right)_t=\frac\psi \theta ,
\ \left(\dsdx \psi q \right)_t=0,\ {\rm with}\ 
\left\{
\begin{array}{l}
q=1,\ \mu=0 \\
q=n,\ \mu=\mu_{\rm b},
\end{array}
\right.
\end{equation}
here $\mu_{\rm b}$ is computed by \eq{eq:var1} with
the value of $M_{\rm b}$ obtained at the bottom of the atmosphere.
The set of equations to be solved on the {\em
equidistant} grid, $q_i\equiv i,\,i=1,\ldots,n$ is therefore:

\begin{equation}\label{eq:syq}\displaystyle
\left\{
\begin{array}{l} \vspace{.1in}
\dsdx\xi q=\left[
-\frac{3G}{8\pi}\left(\frac {M_\odot }{R^2_\odot }\right)^2
\left(\frac{\mu}{\zeta}\right)^2+
\frac{M_\odot}{4\pi R_\odot}\sqrt{ \frac{\mu}{\zeta} } \Omega^2 \right]
\frac{e^{-\xi}\psi}\theta 
\\ \vspace{.1in}
\dsdx\eta q=\dsdx\xi q \nabla \\ \vspace{.1in}
\dsdx\zeta q=\frac 3{4\pi}\frac {M_\odot } {R^3_\odot }
\frac 1\rho \sqrt{\frac{\mu}{\zeta}}\frac\psi\theta \\
\vspace{.1in}
\dsdx\lambda q=\frac{M_\odot }{L_\odot }\sqrt{\frac{\mu}{\lambda}}
\left[\epsilon -\left(\dsdx U t \right)_{\mu_{\rm d}} 
+\frac{\exp{\xi}}{\rho^2}\left(\dsdx 
\rho t \right)_{\mu_{\rm d}}\right]
\frac\psi\theta \\ \vspace{.1in}
\dsdx\mu q=\frac\psi\theta\\ \vspace{.1in}
\dsdx\psi q=0 \\
\end{array}
\right.
\end{equation}
The initial and boundary conditions are straightforwardly
expressed in terms of $q$. Note that the 
derivative, with respect to time, of the specific internal energy and of the density 
should be taken along directions d$M_{\rm d}=0$, \ie, d$\mu_{\rm d}=0$, taking 
into account the change of mass due to the mass defect.
Indeed, the solution, by the damped Newton Raphson scheme, of the non-linear 
\eq{eq:syq}, necessitates the knowledge 
of the derivatives of $\rho$, $\epsilon$ and $U$, with 
respect to the mass, taking into account the changes of chemical composition;
that is the more restrictive conditions imposed by the use of this automatic 
allotment of mesh points; it is also a consequence of the fact that the 
integro-differential problem, 
\eq{eq:eq1}, cannot be solved as a whole since,
along the evolution, convective zones appear, disappear, cede or recede with
not fixed limits.
Moreover with an equidistant mesh in $q$, differences between close numbers
are avoided; that is particularly sensitive in the external part of
the envelope where the changes of mass between two 
adjacent grid points are very small.

\paragraph{Choice of $Q$.}
$Q$ should be a strictly monotonous, two times differentiable, 
function, and as simple as possible. By experiments, 
it has been found\footnote{For models of interiors of giant planets (Guillot \&
Morel, \cite{gm}) another form for $Q$ is used.} that:
\begin{equation}\label{eq:qrp}
 Q=\Delta_\xi\xi+\Delta_\mu\mu
\end{equation}
is, in all, the most convenient form; for the two "distribution factors"
$\Delta_\xi$ and $\Delta_\mu$, the {\em heuristic} values: 
\[
\Delta_\xi\equiv -1,\ \Delta_\mu\equiv 15,
\]
are close to those used by Eggleton (\cite{eg}). The function 
$\theta(\mu,t)$ can now be explicitly calculated using \eq{eq:theta}
and \eq{eq:syq}.
In the core, the pressure gradient is not large and the mesh refinement is 
monitored by the changes of $\mu$, \ie, the mass; while, in the outermost part of 
the envelope the repartition function is controlled by the changes of $\xi$, \ie, 
the pressure, due to its large gradient; there, from a grid point to the next, 
the mass changes are very small, typically $\sim\Delta M\sol 10^{-10} M_\odot$,
even $\sim\Delta M\sol 10^{-12} M_\odot$, while for, the pressure, the changes are of 
the order of $\sim 10\%$ (with $C(t)\sim 0.1$, see next paragraph); a similar
situation occurs on the neighborhood of shell sources.
 
\paragraph{Optimization of the number of grid points.}
Between a PMS initial model and an evolved model at the beginning of the 
$^4$He burning, the central pressure is magnified by more than 
30, that also affects the distribution constant $C(t)$; the accuracy is ensured 
whatever the age is, if the repartition constant, is kept almost fixed with 
respect to time; that is done by increasing (or decreasing) 
the total number of shells in such a way that $C(t)$ remains within $\pm 
2\%$, of its initial value. With $C(t)\simeq 0.1$, $\Delta_\xi=-1$ and 
$\Delta_\mu=15$, the relative change in pressure within a shell is typically 
10\% and the number of zones in a PMS initial solar model is of 
the order of 250, it increases to $\sim 450$ at present solar age and more 
than $\sim 800$ at the onset of the helium flash.

\paragraph{Setting a grid point on a LMR.} The algorithm of automatic 
location of grid points has been extended in order to shift automatically a grid 
point in a 
close vicinity of each LMR; the method is described in \App{sec:set}.

\subsubsection{Boundary conditions at center}\label{sec:bccn}
With respect to the independent variable $q$, the center corresponds to $q=1$,
the three inner boundary conditions are simply written:
\begin{displaymath}
\zeta(1,t)=0,\ \lambda(1,t)=0,\ \mu(1,t)=0.
\end{displaymath}

\subsubsection{External boundary conditions}\label{sec:ebcn}
As seen in \sec{sec:ebc}, at time $t$, the bottom of the atmosphere is connected
at the outermost boundary of the envelope by the three functions: pressure
$P_{\rm b}(R,L,t)$, temperature $T_{\rm b}(R,L,t)$ and mass $M_{\rm b}(R,L,t)$
of the radius, luminosity and time; they are derived from the outer boundary 
conditions through a reconstitution of the atmosphere. 
With the variables of \eq{eq:var1}, at $q=n$, the external 
boundary conditions for the envelope are written:
\begin{equation}\label{eq:ext}
 \begin{array}{l}
  \xi(n,t)-\xi_{\rm b}(\zeta_{\rm b},\lambda_{\rm b},t)=0 \\
  \eta(n,t)-\eta_{\rm b}(\zeta_{\rm b},\lambda_{\rm b},t)=0 \\
  \mu(n,t)-\mu_{\rm b}(\zeta_{\rm b},\lambda_{\rm b},t)=0
 \end{array}
\end{equation}
\noi where $\xi_{\rm b}$, $\eta_{\rm b}$ and $\mu_{\rm b}$ are the 
values, at time $t$ of, respectively, $\xi$, $\eta$ and $\mu$ at the bottom of 
the atmosphere, of radius $R_{\rm b}=\sqrt{\zeta_{\rm b}}R_\odot $ and 
luminosity 
$L_{\rm b}=\sqrt{\lambda_{\rm b}^3}L_\odot$, where the diffusion approximation 
becomes valid.
Hence the grid adjustment solves trivially the problem of the open limit.

\paragraph{Simplified external boundary conditions.}
In that case, $M_{\rm b}\equiv M_\star$ and $\mu_{\rm b}\equiv 1$; $P_{\rm b}$ 
and $T_{\rm b}$, are derived from \eq{eq:sext} with 
$R_\star\equiv R_{\rm b}=R_\odot\sqrt{\zeta(n,t)}$ and 
$L_\star\equiv L_{\rm b}=
L_\odot\sqrt{\lambda^3(n,t)}$. 
\noi The Rosseland mean opacity  $\kappa$ is a function of density and
temperature, therefore the non-linear 
system of \eq{eq:ext} must be solved by iterations for $\xi_{\rm b}$,
$\eta_{\rm b}$, $\mu_{\rm b}$ and for their derivatives with respect to $\zeta$
and $\lambda$.

\paragraph{Precise external boundary conditions.}
As described in \Sec{sec:ebc}, the precise reconstitution of an atmosphere consists
in a differential problem, with boundary conditions at three different levels:
$\tau_b\goa 10$, $\tau_{\rm min}\loa 10^{-4}$ and $0.3\loa\tau_\star\loa 0.6$; 
recall that, at $\tau=\tau_\star(t)$, defined by \eq{eq:star1}, there is a {\em 
open} inner limit; the solution given to this numerical challenge is described 
in \App{sec:pbc}.


\subsubsection{Gravothermal energy}\label{sec:tds}
The gravothermal energy is written\footnote{As pointed out by Strittmatter 
\Etal{st}, Reiter \Etal{re1},
this formulation is correct while, in Eq. 4.27 of Kippenhahn \& Weigert 
(\cite{kw}), the changes of chemicals are ignored.}: 
\begin{equation}\label{eq:tds}
 -\epsilon_{\rm G}=\dsdx Ut-\frac P{\rho^2}\dsdx\rho t.
\end{equation}
For stability  purposes, at each collocation point, $\epsilon_{\rm G}$ is
discretized by the backward difference formula of first order:
\begin{equation}\label{eq:tdsd}
 -\epsilon_{\rm G}=\frac{U^{t+\Delta t}-U^t}{\Delta t}-\frac  
 P{\rho^2}\frac{\rho^{t+\Delta t}-\rho^t}{\Delta t}+O(\Delta t^2).
\end{equation}
The use of this low order formula is justified by the fact
that the entropy variations are negligible on the main-sequence, otherwise 
the models are so complicated that
high order schemes can lead to instabilities and therefore are not recommended. 
In the plane ($M\times t$), Eq.~(\ref{eq:tds}) must be written along
directions of constant mass, \ie, including the mass defects,
d$M_{\rm d}=0$ or d$\mu_{\rm d}=0$; at 
time $t$, density and specific internal energy are derived, via 
EOS, from pressure, temperature and chemical 
composition (see Sect.~\ref{sec:chim}) obtained through the {\em non-linear} 
inverse interpolation scheme ``$\hookrightarrow$'':
\[
 M_{\rm d}\rightarrow\mu_{\rm d}\rightarrow\left\{
 \begin{array}{ll}
  \hookrightarrow q & \mbox{$\rightarrow(\xi,\eta) \rightarrow(P,T) $,} \\
  \rightarrow\nu    & \mbox{$\rightarrow{\cal X}$,}
 \end{array}
 \right.
\]
\noi with $\nu\equiv (M/M_\odot )^{2/3}$ (see \sec{sec:chim}).

\subsubsection{An eulerian set of variables}\label{sec:eul}
Note that, with $(R,t)$ as an eulerian independent variable, \eq{eq:eq1} have
no singularity but, as the stellar radius varies with respect to time, the 
external limit, at 
$M=M_{\rm b}$, is an open limit; though that is not a difficulty for the 
numerical 
method of integration used -- see \sec{sec:aloc} --, the lagrangian set of 
variables 
\eq{eq:var1}, is preferred owing to its highest numerical accuracy as 
demonstrated \App{sec:ap1}. With the 
exponent 2/3 for $\lambda$ , in \eq{eq:var1}, the luminosity must remain 
non-negative, \ie, $\lambda\geq 0$; otherwise, though less accurate, the 
following 
set of eulerian variables ought to be employed:
\begin{equation}\label{eq:eq5}
 \begin{array}{l}
\xi=\ln P,\ \eta=\ln T,\ \mu=\frac M{M_\odot },\ \zeta=\frac R{R_\odot },
\ \lambda=\frac L{L_\odot }.
 \end{array}
\end{equation}
\medskip
Similarly, for the eulerian variables, one defines:
$\psi(t)=\left(\dsdx Q q\right)_\zeta$ and
$\theta (\zeta,t)=\left(\dsdx Q{\zeta}\right)_t$.
\noi Equations, initial and boundary conditions are straightforwardly derived.

\subsection{Numerical solution of the "evolution of chemicals"}\label{sec:evol}
To be consistent it is desirable to be able to solve, with algorithms of same 
orders, the two differential problems, \ie, the "quasi-static equilibrium" and the 
"evolution of chemicals".
One needs also a suitable monitoring of the time step derived, if possible,
from an estimate of the accuracy of the integration.

\subsubsection{Interpolation of chemicals}\label{sec:chim}
Due to the change of mass caused by the mass defect and mass loss, the
chemicals need to be  known at any point in the star. Here, 
with the spline-collocation method employed for the numerical integration of
the quasi-static equilibrium, the collocation points, where the 
discretized differential equation are fulfilled, differ from the grid points; 
moreover, due to the grid adjustment, the collocation points have not fixed 
values in mass; therefore an interpolation algorithm is needed for the 
knowledge of the chemical composition at any point of the interval 
$[0,M_\star(t)]$. In fact, interpolation is needed as far as grid points are 
moved
or added during the evolution; on the other hand, 
interpolations create numerical diffusion which, in some extent, 
stabilizes the numerical solution. The interpolation scheme needs to be, at
least, as precise as the integration algorithm.
Here one seeks the distribution of chemicals as piecewise polynomials
on a lagrangian mesh. Due to the assumption of 
spherical symmetry, in the vicinity of the center the gradients of 
abundances are written:
\[
 \dsdx{X_i}R=O(R),\ i=1,\ldots,n_{\rm X}.
\]
\noi Therefore
\[
 \dsdx{X_i}M=\dsdx{X_i}R\dsdx RM=O(R^{-1}),
\]
\noi with respect to the mass, the first derivative of the abundances
have a singularity
at the center so it cannot be interpolated accurately with respect to $M$
using the variable such that $\nu\equiv(M/M_\odot)^{2/3}$:
\[
 \dsdx{X_i}{\nu}=\dsdx{X_i}R\dsdx R{\nu}=\frac{O(R)}{O(R)}=O(1),
\]
\noi the singularities are avoided.
With diffusion (see \sec{sec:mixd}), 
the solution of the diffusion equation agrees with the piecewise
polynomial of interpolation;
otherwise piecewise polynomial interpolation, with order $h_c$, $2\leq h_c\leq 
4$, with respect to $\nu$ is used 
(see Appendix~\ref{sec:int}). An {\sl ad hoc} choice of the B-spline basis
allows to take the discontinuities into account. 
Owing to its linearity, the interpolation
scheme is conservative for any linear set of the unknowns -- a proof is given 
in Appendix~\ref{sec:cpw}; therefore
the total charge, total baryon number and baryon numbers involved in the PP
and CNO cycles remain constant, within roundoff errors.

Due to (i) the splitting in two steps of the integration of the whole problem, 
and (ii) to the ability of restoring the solution at any point, a grid can be 
especially designed for the chemicals; this net is refined in the inner parts 
(recall that thermonuclear reaction rates have 
a high power law dependence with respect to the temperature)
where the nuclear reactions are active.

\paragraph{Management of discontinuities.}
As noticed previously, a long-lasting discontinuity on chemical composition is
an unphysical situation; therefore, when the diffusion is ignored, if no 
physical process is explicitly introduced, the discontinuities are smoothed by 
the 
numeric, the lower the order of the numerical scheme is, the more efficient is 
the smoothing. Satisfactory results have been obtained
using, for a given time step, the mass points designed for the quasi-static
problem, and for the next time step, mass points located at half distance
between two neighbouring mass points designed for the quasi-static
problem; ensuring, however, that the 
mass step, for the chemical composition is, at least, greater than
$\Delta m \geq 5\,10^{-4}M_\odot$, except for each CZ where a minimum of 10 
points is required.
As far as a convective core increases, the chemicals undergo discontinuities at
its limit, while, when it recedes, at any point localized in the zone between 
the
previous and the new limits of the core, the values of abundances not only 
depend on the local temperature and density but, also, on how long that
point has lasted in the mixed receding core; a similar situation occurs as
soon as a 
MZ recedes from a radiative zone. Moreover in radiative zones, when the
diffusion is 
ignored, though violating the physics, the 
discontinuities in chemicals {\em formally} stay so far they are not, again, 
embedded in a new extent of a convective zone. It 
is difficult to mimic in details and precisely all these tricky processes.
In CESAM, 
the abundances at any point localized between the limit of a core and its
location at the former time step, are obtained by linear 
interpolation, with respect to $\nu$, between their values at the new and at the 
previous locations of 
the limit.

\subsubsection{Standard evolution}
\label{sec:mixc}
In the radiative zones, without diffusion, the equations to 
be solved are written:
\begin{equation}\label{eq:chim}
 \dsdx {X_{\rm i}}t =  \Psi_i,\ 1\leq i\leq n_{\rm X}.
\end{equation}
In MZ, the chemical composition is homogeneous and the equations for the mean
abundances can be written (see sect.~\ref{sec:mixisd}):
\begin{eqnarray}\label{eq:dx2}
 \ddx {\bar{X_i}}t&=&\int_{\rm MZ}\Psi_i(P,T,\bar{\cal X},t)
 {\rm d} M \bigg/  \int_{\rm MZ} {\rm d} M \nonumber \\
 &=& \int_{\rm MZ}\sqrt\nu\Psi_i(P,T,\bar{\cal X},t)
 {\rm d} \nu \bigg/  \int_{\rm MZ} \sqrt{\nu}{\rm d} \nu,
\end{eqnarray}
therefore it is an integro-differential problem.
Another numerical difficulty results from large ratios between the
characteristic evolutionary time scales involved;
they can be estimated by the ratio between the eigenvalues, $\lambda_{\rm min}$
and $\lambda_{\rm max}$, of minimum and maximum norms of the jacobian matrix:
\[
J=\left(\dsdx{\Psi_i}{X_j}\right)_{i,j=1,\ldots,n_{\rm X}}
\]
\noi Typically, with the physical conditions at
the center of the present Sun, $T\sim 15\, 10^6$K, $\rho\sim 150$g cm$^{-3}$, 
one has:
\begin{displaymath}
\frac{|\lambda_{\rm max}|}{|\lambda_{\rm min}|} \ga 10^{18}.
\end{displaymath}
Such differential problems are called "stiff"
(Gear, \cite{ge}; Hairer \& Wanner, \cite{hw}); special algorithms have been 
developed for
their numerical solution; they ensure numerical stability even if the
time step is larger than the smallest characteristic time scale. However,
it is not possible to have accurate solution
for all the variables, regardless of the size of the time step;
therefore, owing to the stability of the scheme, the numerical
errors are damped out, a given accuracy being ensured for variables of interest.
Furthermore, chemical abundances being positive numbers, oscillations around
zero (as observed with the trapezoidal rule) are unsatisfactory.
The so-called "L-stable" schemes (Hairer \& Wanner, \loc)
have good stability properties without oscillations; they are suitable
for the integration of the evolution of chemicals abundances.
Equations~(\ref{eq:chim}) 
and~(\ref{eq:dx2}) and their numerical counterparts are non-linear; they 
are also implicit, as the L-stable schemes are.
The L-stable Implicit Runge Kutta (IRK) Lobatto IIIC 
formula with orders\footnote{according to standard definitions the local
error is $O^{p+1}$.} $p=$1, 2 and 4 are available in CESAM; 
their coefficients are reproduced table~\ref{tab:lob} (\App{sec:sirk}), $p=1$ is the standard 
Euler's backward scheme (Hairer \& Wanner, \loc).
For the Lobatto IIIC formulas with order $p$ greater than two, values for 
the 
temperature and density are needed at intermediate time levels. They are
estimated by interpolations of order four from successive models.

\paragraph{Stable integration with IRK formulas.} Emphasis is made on the 
fact that, with stiff problems it is of great importance to use special 
algorithms specially designed for IRK formula; those used in CESAM are described 
in \App{sec:sirk}.


\paragraph{Numerical control of the accuracy.}
The comparison of solutions given by two IRK formula
which differ by one order of accuracy, \ie, the so-called Fehlberg method (Stoer
\& Bulirsch, 
\cite{sb}), allows an estimate $\varepsilon_e$ of the numerical accuracy.
Here the IRK formulas Radau IIA (Hairer \& Wanner, \loc, \S IV.8)
are used in connexion with Lobatto IIIC formulas.
Let $\varepsilon_r$ be the value required for the relative precision; an
optimal value $\Delta t_{\rm opt}$ for the next time step is written:
\begin{displaymath}
 \Delta t_{\rm opt}=\left(.9\Delta t\frac{\varepsilon_r}{\varepsilon_e}
 \right)^{\frac 1{1+p}},
\end{displaymath}
here $p$ is the order of the IRK formula.
It has been observed that the robustness of the scheme is improved if
only small changes for the time step are allowed for; therefore the estimate of
the new time step $\Delta t_{\rm new}$ is taken as:
\begin{equation}\label{eq:ts}
 \Delta t_{\rm new}=\max\left( 0.8\Delta t,\min\left[1.2\Delta t,
 \Delta t_{\rm opt} \right] \right)
\end{equation}

This precise control of the numerical accuracy, practically, doubles the 
computing time; it is prohibitive in most cases, therefore the time step is 
simply adjusted in such a way that the 
changes of the abundances remain within fixed limits
$\Delta_i,\ i=1,\ldots,n_X$.

\subsubsection{Evolution with diffusion}
\label{sec:mixd}
With diffusion, for every $\nu\in [0,\nu_{\rm b}]$ 
($\nu_{\rm b}\equiv(M_{\rm b}/M_\odot)^{2/3}$), the set of equations to be 
solved is written:
\begin{equation}\label{eq:diff}
 \dsdx {X_{\rm i}}t =-\dsdx{F_i}M+\Psi_i=-\frac2{3M_\odot 
\sqrt\nu}\dsdx{F_i}\nu
 +\Psi_i,
\end{equation}
for $1\leq i\leq n_{\rm X}$. 
As seen \sec{sec:mixiad}, the mixing in MZ is made by turbulent
diffusion. Since \eq{eq:diff} holds everywhere, the evolution 
of the chemical composition is no longer an integro-differential problem but a 
differential problem with
boundary conditions given by \eq{eq:bce} and \eq{eq:bcn}; it is a mixed 
parabolic/hyperbolic problem.
At each LMR, the abundances $X_i$ and the fluxes $F_i$ are continuous
functions 
with discontinuous first derivatives owing to the jumps of $d_{\rm M}$ (see 
\sec{sec:mixiad}).
The method of integration of the diffusion equation, written in finite-elements 
form, is described in \App{sec:ided}.

\paragraph{Efficiency of the mixing by diffusion.}
\Fig{fig:d_vson} plots, with respect to the mass fraction, the 
relative differences in
sound velocity for two calibrated solar models calculated,
respectively, with standard mixing and mixing by 
diffusion  ($D^\star_{i,j}+D_{\rm T}=0$, $D_{\rm M}=10^{13}$ and $v_i=0$,
$\forall i,j$).
The differences, at the level of $10^{-5}$, show the similar efficiency of
both approaches.
\Fig{fig:abon} plots, with respect to the radius fraction,
the normalized abundances of $^2$D, $^3$He, $^7$Li, $^7$Be, $^{12}$C at the 
onset of
the main sequence for a non-standard solar model, evolved from PMS,
including  microscopic diffusion according to Michaud \& Proffitt (\cite{mi});
at that age (45My) the convective core of the young Sun recedes, while the CZ
is close to its present day location. One emphasizes on the well marked
drop of the gradient of $^2$D and on the smooth profile of $^7$Li at 
the limit of the
CZ; due to the mixing, in the
convective core, the elements have a constant abundances except $^8$Be
since its nuclear time is of the order of the mixing time ($\sim$100 days). As 
seen, the
Petrov-Galerkin's solution is stable even with strong jumps of more than 
thirteen decades in the diffusion coefficients profiles.
\begin{figure}
\defbox{\psfig{figure=diff_vson.ps,%
height=7cm,bbllx=0.5cm,bblly=0.5cm,bburx=27.5cm,bbury=20.5cm,clip=}}
\caption{Relative difference of the sound velocity
between two calibrated solar models
calculated with the two kinds of mixing: standard instantaneous mixing and 
diffusion. The low level of the differences shows that the two kinds of mixing 
have similar effects. The wiggly behavior for radius lesser than 
$0.1M_{\rm tot}$ is a fossil signature of the discontinuities of
chemicals at the limit of the convective core of the young sun; they are
almost not smoothed because these calculations employed a low 
dissipative scheme for the interpolation of chemicals. 
}
\label{fig:d_vson}
\end{figure}


\begin{figure}
\defbox{\psfig{figure=des_abon.ps,%
height=7cm,bbllx=0.5cm,bblly=0.5cm,bburx=27.5cm,bbury=20.5cm,clip=}}
\caption{Normalized abundances in a $1M_\odot$ non-standard  model at the end
of the PMS. At that age (45My) the convective core of the young Sun recedes and
the CZ is closed to its present day location. At LMR the solution is stable
even if
jumps of more than thirteen decades affect the diffusion coefficients. The 
maxima respectively are:
$^2$D=$1.6\,10^{-17}$, $^3$He=$3.1\,10^{-4}$, $^7$Li=$4.1\,10^{-9}$,
$^7$Be=$1.7\,10^{-12}$, $^{12}$C=$3.1\,10^{-3}$.
}
\label{fig:abon}
\end{figure}

\section{General features for the implementation of the code}\label{sec:gf}
The flow chart of CESAM takes advantage of two spaces: a functional 
space (B-spline), where the differential equations are integrated with the 
mathematical formalism and, a physical space, where the 
equations are written regarless of the method employed for their solution; it 
ensues a modular structure which allows to exploit, with the same algorithm,
several sets of physical data, \eg 
EOS, opacity, thermonuclear reactions, diffusion, atmosphere, \etc
and a few numerical parameters allow to adjust the order and the accuracy of 
the solution.

\subsection{Initializing and updating the solution}\label{sec:upd}
The iterative process is initialized with a
model which is, either, for $t=0$, a starting model of PMS or ZAMS or,
for $t>0$, a previously evolved model;
a new initial model is then computed with the requested open parameters for 
physics, \ie, mass, mixing length, \etc, and for numeric, \ie, number of shells, 
order of piecewise polynomials \etc
Along the evolution, the calculation of a new model is initialized with the 
previous model without extrapolation.
When the gravothermal energy is the leading term in the energy equation, such
as occurs during the PMS phase or with shell sources,
the values for density and specific internal energy need to be accurately 
initialized (Kippenhahn \etal{ks}); since the approximate model is the previous
one, the initial value for the gravothermal energy is zero, the energy 
equation 
is strongly in error and the Newton-Raphson scheme does not converge. The
remedy of H\" arm \& Schwarzchild (\cite{hs}), consists in taking the
approximate model as the previous one shifted in mass; in the B-splines space, 
this technique is not
easy to be worked out and, by experience, it does not work for PMS. A trivial 
solution has been found: for the first Newton-Raphson 
iterate,
{\em but only for this first}, the gravothermal energy is taken equal to its 
value in the former model; hence the energy equation remains in almost closed 
form (recall that 
is one of the necessary conditions for the convergence of the Newton-Raphson 
method).

\subsection{Time step controls}\label{sec:tsc}
There are three levels of time step controls; as soon as one of the
following criteria is not fulfilled the time step is divided by two,
otherwise, for the next model it increases, at most, by 20\%:
\begin{enumerate}
 \item inadequate accuracy of the solution or no convergence
 of the IRK scheme, after more than 12 Newton-Raphson iterations,
 \item slow convergence or divergence of the Newton-Raphson scheme
  used for the two points boundary problem,
  \item too large changes in chemicals,
 \item too large gravothermal energy change.
\end{enumerate}



\subsection{Sets of numerical parameters}\label{sec:tv}
The amount of calculations increases with the accuracy to be achieved, for some 
applications precise models are not necessarily needed, \eg with uncertain or 
unstable physics; for others applications, high accuracy is reachable and is of 
necessity, \eg calibrated solar models used for 
the calculations of frequencies; therefore the 
numerical parameters have to be adjusted to the aims towards the model is 
directed. \tab{tab:np} gathers three sets of numerical parameters of common use: 
i) (la) {\em low accuracy}, for tests, uncertain or unstable physics, ii) (sa) 
{\em standard accuracy}, iii) (Sa) {\em Solar accuracy}, for the 
calculation of precise solar models. 

\begin{table}
\caption[ ]{
Sets of numerical parameters employed for the calculation of models with~:
{\em low accuracy} (la), {\em standard accuracy} (sa), {\em Solar accuracy} (Sa) 
(the other notations are defined \sec{sec:tv}). 
}
\begin{tabular}{lllll} \hline
                   & (la) & (sa) & (Sa) &\\ 
\hline
$C$                & 0.11 & 0.1  & 0.06 &\\

$h_c$              & 2(2)       & 3(3)       & 4(3)       &\\ 
$p$                & 1          & 2          & 2          &\\
$(\Delta X)_{\max}$& 0.3        & 0.2        & 0.1        &\\
$(\Delta t)_{\max}$& 200        & 100        & 50         &\\

$h_q$              & 2(2)       & 3(4)       & 3(4)       &\\
$\Delta_{\rm N}$   &$5\,10^{-3}$&$1\,10^{-4}$&$5\,10^{-6}$&\\
$\Delta_g$         & 1.0        & 0.5        & 0.1        &\\

$n_{\rm a}$        & 25         & 30         & 50         &\\
$n_\star$          & 15         & 20         &  30        &\\ \\

$n_{\rm PMS}$      & 250        & 270        &  450       &\\
$n_{\rm ZAMS}$     & 390        & 420        &  700       &\\
$n_\odot$          & 395        & 430        & 720        &\\
\hline

 \end{tabular}\label{tab:np}
\end{table}
\noi  In \tab{tab:np}, $C$ is the repartition constant, $h_c$ the order of 
piecewise polynomials employed for interpolation of chemicals with (in 
parenthesis) their values with diffusion; $p$ is the order of the IRK scheme; 
$(\Delta X)_{\max}$ is the maximum relative change allowed for $^1$H and $^4$He 
during a time step of maximum length $(\Delta t)_{\max}$; $h_q$ is the order of 
piecewise polynomials of the solution with (in parenthesis), 
their orders magnified by super 
convergence (see \sec{sec:scv}); $\Delta_{\rm N}$ is the closure used for the damped Newton-Raphson 
scheme; $\Delta_g$ is the maximum relative change allowed to the gravothermal 
energy in the core during a time step; $n_{\rm a}$ and $n_\star$ respectively 
are the total number of shells in the atmosphere and between the bottom of the 
atmosphere and the level where the radius of the star is located; $n_{\rm PMS}$ 
$n_{\rm ZAMS}$, $n_\odot$ respectively are typical numbers of shells in solar 
model at the onset of the PMS, of the ZAMS and for the present solar age. 

\subsection{Calibration of solar models}\label{sec:cal}
The solar models are calibrated by 
adjusting the ratio $\alpha\equiv l/H_{\rm p}$ of the mixing-length to the 
pressure scale height, 
the initial mass fraction $X$ of 
hydrogen and the initial mass fraction $Z/X$ of heavy elements to hydrogen in 
order to have, at present day, the luminosity $L_\odot$, the radius $R_\odot$ 
and the mass fraction ${\cal R}_{Z\odot}\equiv(Z/X)_{t=t_\odot}$. According to 
Guenter \Etal{gu} the most 
recent determinations for the solar luminosity, radius, mass and age are:
$L_\odot=3.8515\,10^{33}$erg sec$^{-1}$, 
$R_\odot=6.9598\,10^{10}$cm, 
$M_\odot=1.9891\,10^{33}$g, and $t_\odot=4.52$Gy.
The value of the ratio of heavy elements to hydrogen, at present day is, 
${\cal R}_{Z\odot}=0.0245$ (Grevesse \& Noels, \cite{gr}; Noels \etal{no}).

With the microscopic diffusion coefficients of Michaud \& Proffitt (\cite{mi}), the 
following dependences on the initial PMS values of the unknown parameters has 
been 
found:
\[
\begin{array}{l}
{\left( \begin{array}{lll}
-0.6471 & -20.18 & -98.51 \\ 
0.07193 & -42.50 & -222.9 \\ 
2.758\,10^{-4} & 1.476\,10^{-2} & 1.016
\end{array}\right)
\left( \begin{array}{l}
\Delta\alpha \\ 
\Delta X \\ 
\Delta {\cal R}_{Z0}
\end{array}\right) =} \\ \\
\qquad\qquad{=\left( \begin{array}{l}
\Delta R\,10^{10} \\ 
\Delta L\,10^{33} \\ 
\Delta {\cal R}_{Z\odot}
\end{array}\right)}\end{array}
\]
\noi with ${\cal R}_{Z0}$ as the value of $Z/X$ at $t=0$ and ${\cal 
R}_{Z\odot}$ its value at present solar age.
Along the calibration process the matrix of the dependence can be updated, using 
the one rank Broyden's method (Conte \& de\,Boor, \cite{cdb}, p. 222)

 



\section{Estimates of the numerical accuracy}\label{sec:test}
For results concerning external numerical tests using the simplified physics
of the Solar Model Comparison project (SMC)
 see Christensen-Dalsgaard, (\cite{cd1}, \cite{cd2}).

\subsection{Physics used for the numerical tests}\label{sec:ctes}
Otherwise stated, the following physics is employed in numerical tests:

\begin{itemize}
\item EOS: EFF formalism (Eggleton \etal{eff}).

\item Opacities: 4D linear $(T_6,A,X,Z)$ interpolation of OPAL radiative 
opacity 
data (Rogers \& Iglesias, \cite{ri}) extended with Kurucz's (\cite{kz}) opacity 
data for low temperatures; as usual, with OPAL opacity data, $T_6\equiv 
T\times10^{-6}, A\equiv \ln\rho/T_6^3$.

\item Convection: at each level, the Schwarzchild's criterion is employed to 
decide if the energy transport is made by radiation or by convection. The basic
formulation of the classical mixing length theory according to Henyey \Etal{hv}
is used for the calculation of the convective gradient; it also 
takes into account the optical thickness of the convective elements.

\item Thermonuclear reactions: the nuclear network contains the following twelve
species which enter into the most important nuclear reactions of the PP, CNO
and the  beginning of  the  $3\alpha$ cycles chain (Clayton,  \cite{cy}): $^1$H,
$^2$D, $^3$He, $^4$He, $^7$Li, $^8$Be, $^{12}$C, $^{13}$C, 
$^{14}$N, $^{15}$N, $^{16}$O, $^{17}$O; in some models, the species $^2$D, $^7$Li, $^8$Be are 
assumed to be at equilibrium; the relevant nuclear reactions are interpolated 
from the tabulated rates of Caughlan and Fowler (\cite{cf}); weak screening is 
assumed; the initial abundances are taken from Anders \& Grevesse~(\cite{ag}) 
meteoric abundances. 

\item Atmosphere: single  shell approximation  or  restoration of an  atmosphere
using a ${\cal T}(\tau)$ law derived from Kurucz's model atmospheres (Morel \etal{mp4}). 

\item Diffusion coefficient: if diffusion is  not ignored,  the  microscopic
diffusion coefficients of Michaud \& Proffitt (\cite{mi}) are employed. 

\item Angular momentum: the rotation is ignored.

\item Mass loss: the mass loss is assumed to be zero and the mass defect is 
ignored.

\end{itemize}

\subsection{Superconvergence}\label{sec:scv}
The expansion of the unknown functions on a basis allows to recover 
exactly the numerical solution at any location, \ie, {\em not only at the grid 
points}. Moreover with the de\, Boor's choice of the collocation points (de\, 
Boor,~\cite{db}),
the order of the accuracy of the numerical solution, at the grid points, is 
improved\footnote{only where the true solution is sufficiently smooth.} by {\em 
superconvergence} (see \app{sec:ide}). 
This nice numerical property avoids the Richarson extrapolation 
(Christensen-Dalsgaard, \cite{cd0}) which is, with stiff problems, sometimes 
unstable (Hairer \& Wanner, \cite{hw}). 
That improvement is illustrated \Fig{fig:super}: for a 
calibrated standard solar model, evolved from the PMS to the present age with 
"standard accuracy" (sa), but with a {\em constant number of shells} fixed to 
300, the 
relative accuracy on the pressure gradient $\dsdx PM$:
\[ \left( \dsdx PM +\frac {GM}{4\pi R^4}\right) \bigl/ {\frac{GM}{4\pi R^4}} \]
\noi is plotted for points randomly distributed on the interval $[0,R_\odot]$;
as seen, the relative precision, taking advantage of the super convergence, is 
improved by one order of magnitude.
\fig{fig:super} illustrates also the stability of the numerical solution near 
the center, the possibility of recovering, with full accuracy, the 
numerical solution everywhere and, the significant 
improvement of the accuracy in the envelope due to the superconvergence.
\begin{figure}
\defbox{\psfig{figure=super.ps,%
height=8cm,bbllx=.5cm,bblly=.4cm,bburx=14.5cm,bbury=13.5cm,clip=}}
\caption{Comparison of the relative accuracy on the pressure gradient 
between the numerical solution and its improvement taking advantage of 
the super convergence. Clearly one order of magnitude is acquired. One 
also emphasizes on the good behavior of the solution in the vicinity of the 
center.
}\label{fig:super}
\end{figure}

\subsection{Estimate of the internal accuracy}\label{sec:iacc}
For more than 3500 points randomly distributed
in a model of 1 $M_\odot$, evolved from 
ZAMS to 4550\,My, with "Solar accuracy" (Sa), the relative differences
between left and right hand sides of the four first differential equations of
\eq{eq:eq1}, namely the derivatives $\dsdx PM$, $\dsdx TM$, 
$\dsdx RM$ and $\dsdx LM$, have been computed from the improved solution of
fourth order; the results, plotted \Fig{fig:preci}, allow an estimate of the
internal accuracy achieved. That model includes an atmosphere reconstructed (see 
\sec{sec:ebcn}) and smooth opacity data have been
provided using the bivariate rational spline interpolation scheme of
Houdek \& Rogl (\cite{ho}).

As seen, the differential equations are fulfilled within a relative 
precision better than $\sim 5\,10^{-6}$ except in the neighbourhood of
the LMR at $R\approx 0.72R_{\rm tot}$, and in the atmosphere for
$R\ga 0.9992R_{\rm tot}$. 
At LMR the local drop of the precision is a consequence of the small
discontinuity of the density, caused by the mixing in the convective zone
where the nuclear reactions, though weak, are active; there $\dsdx RM$ is
discontinuous, the solution is only of class ${\cal C}^0$, the fourth order
accuracy cannot be reached; the effect is more sensitive for $\dsdx RM$ 
than for the other gradients, \fig{fig:preci} panel (c). The cause of
the low accuracy $\sim 2\,10^{-3}$ in the restored atmosphere, panel (d),
is a consequence of the non-linear dependence of the {\em natural} variable 
$\tau$, used for the restoration of the atmosphere, with respect to $R$ through 
the ${\cal T}(\tau)$ law and the opacity.

In the model, the closest grid 
point from the center is located at $R\sim1.2\,10^{-2} R_{\rm tot}$;
in the plots, deliberately, the first point was located
at $R\sim3.8\,10^{-4} R_{\rm tot}$, \ie, quite near to the central singularities;
the stability of the numerical solution there is clearly illustrated 
\fig{fig:preci}, panels (a), (b) (c) and (e).

As demonstrated \app{sec:cirk} and \app{sec:cpw}, the algorithm employed for the 
calculation of the chemical composition are conservative; that is illustrated 
\fig{fig:preci} panel (f), where the local departure from the mean value of the 
number of baryons,
$\sum_{i=1}^{n_X} \frac {n_i}{A_i} X_i$
is plotted; $n_i$ and $A_i$ are, respectively, the baryonic number and the 
atomic weight of $X_i$. As seen, the conservation of the number of baryon is 
ensured within the machine accuracy except at LMR, due to the mixing which is 
not made with so high numerical accuracy (see \sec{sec:mixc}) and, also to the 
displacement of the LMR.

\begin{figure*}
\defbox{\psfig{figure=preci_sa.ps,%
height=12cm,bbllx=-1.cm,bblly=-.4cm,bburx=16.5cm,bbury=18.5cm,clip=}}
\caption{Panels (a), (b), (c), (e): relative accuracy on the closure for, 
respectively, the first four differential equations of stellar evolution; panel 
(d): enlargement for outer layers for the pressure gradient; panel (f): relative 
difference on number of baryons. 
These plots are computed, using the improved solution, for more than 3500 
points, randomly distributed in the model; the radius of the innermost point is 
$R\sim3.8\,10^{-4} R_{\rm tot}$.
The accuracy of the solution is better than $\sim 5\,10^{-6}$ except in the 
neighbourhood of the LMR and in the atmosphere.
} 
\label{fig:preci}
\end{figure*}

\subsection{Discontinuities.}\label{sec:dcs}
\fig{fig:vais} illustrates the ability to restore the solution at any location, 
even with discontinuous functions. Due to the convective mixing, in a stellar 
model taking diffusion into account the gradient of the mean molecular weight 
(here noted $\mu$) and the quantity~:
\[
{\cal N}^2\equiv \frac{GM}{R^2}\left(\frac 1{\Gamma_1}\dsdx{\ln P}R-
\dsdx{\ln\rho}R\right), 
\]
are discontinuous at LMR --  $\cal N$, difference between two gradients,
is equivalent to a second derivative -- $\Gamma_1$ is the adiabatic exponent. In 
the neighbourhood of the LMR of a model\footnote{This calibrated solar model was 
evolved, using standard accuracy (sa), from PMS to present solar age with no 
chemical at equilibrium, a reconstructed atmosphere and an under-shooting of 0.2 
pressure scale height.} with microscopic diffusion (Michaud \& Proffitt, 
\cite{mi}), the above discontinuous functions have been derived from the 
solution, taking advantage of the superconvergence; on \fig{fig:vais} more 
than 20 point are inserted between two adjacent grid points. 
Emphasize is made on the following facts~: i) the discontinuity is well marked, 
ii) in its vicinity, even with jumps, 
for the diffusion coefficients, as large as thirteen magnitudes, the solution has
no oscillation, iii) $\cal N$, a second 
derivative, is smooth, iv) a grid point is placed right on the LMR.

\begin{figure}
\defbox{\psfig{figure=des_vais.ps,%
height=7cm,bbllx=-1.5cm,bblly=0.cm,bburx=28cm,bbury=19.5cm,clip=}}
\caption{Plots, with respect to the normalized radius, of 
discontinuous quantities $\cal N$ (dots) and $\partial {\mu^{-1}}/\partial R$ 
(full) -- $\mu$ is the mean molecular weight --  in the neighbourhood of the LMR
of a solar model with microscopic diffusion; on the $\cal N$ curve,
between two grid points $(\bullet)$, each dot corresponds to a {\em 
calculated} and not to an interpolated value of $\cal N$. 
} 
\label{fig:vais}
\end{figure}


\subsection{Comparison of evolutionary tracks}\label{sec:cev}
Two models of $4M_\odot$ have been evolved starting from PMS to the onset of
the $^4$He burning phase; these two models differ only by the numerical
accuracy achieved. The first one, is computed with low accuracy (la), the
second with 
standard accuracy (sa), see \tab{tab:np}; the numbers 
of, time steps, shells at the onsets of PMS, ZAMS, POST and $3\alpha$ cycle
burning phase, are given \tab{tab:tab1}. The
evolutionary tracks are notably stable, even for the rugged solution obtained
with (la); notice also that, on the main sequence,
the two tracks are superimposed.

\begin{table}
\caption[]{Typical number of shells used in the calculations of 
$4M_\odot$ models 
with, respectively, low (la) and standard (sa) accuracy.}  
\label{tab:tab1}
\begin{tabular}{lll} \hline
& (la) & (sa) \\ \hline
$n_{\Delta t}$ & 103 & 137 \\
$n_{\rm PMS}$ & 250 & 270\\
$n_{\rm ZAMS}$ & 420 & 460 \\
$n_{\rm POST}$ & 430 & 480\\
$n_{3\alpha}$ & 480 & 530\\ \hline
\end{tabular}
\end{table}

\begin{figure*}
\defbox{\psfig{figure=hr_pub.ps,%
height=9cm,bbllx=0.cm,bblly=0.cm,bburx=28.cm,bbury=19.5cm,clip=}}
\caption{Evolutionary tracks of two $4M_\odot$ models computed with low
(dotted) and super (full) accuracy from PMS initial conditions to the onset of
the $3\alpha$ cycle. Large time steps are responsible of the angular behavior
of the solution of low accuracy; the two tracks do not significantly
differ, they are superimposed on the main sequence.
}\label{fig:hr}
\end{figure*}

\section{Discussion}\label{sec:dis}
The principal novel feature of the CESAM code, devoted to the calculation of 
stellar evolution, is the assumption that the dependent 
variables can be approached by piecewise polynomials which are, in turn, 
projected on a convenient basis for the numerical analysis; that allows
to handle the 
discontinuities without too much difficulties, to have a scheme with an order
of accuracy having only one parameter dependence and  the 
possibility of improving the numerical solution via the  super convergence.

Other unusual features are the use of i) an automatic location of the numerical
mesh with a variable total number of points and adjustment of grid points
in close neighborhoods of each limit of convective zones, ii) a stiffly stable 
implicit Runge-Kutta scheme with orders up to four for the computation of
the chemicals change, iii) finite elements approach for the calculation 
of diffusion, using Petrov-Galerkin formalism with a large turbulent diffusion 
for mixing the chemicals in convective zones, iv) precise reconstitution of the 
atmosphere with the stellar radius located at the level where the 
local temperature is equal to the effective temperature, v) mass changes due to
the mass loss and mass defect are also allowed for.

The linearized equations are solved using the 
damped Newton-Raphson that allows to follow evolution with good time steps
and reduced nets, even with shell sources. Although the set of lagrangian 
variables used so far allows to reach the best numerical accuracy, it does not 
permit to pursue the evolution towards more advanced phases if negative values 
of the local luminosity are expected; in such cases the calculations are 
possible with a numerically, less accurate set of eulerian variables.

The most constraining features are i) the use of the pressure, instead of 
the density, as state variable and ii) the necessity, for all the physical 
quantities, of numerically consistent derivatives with respect to the variables
of state, \ie, pressure, temperature and chemical composition.

The extension of this code towards more advanced stages of stellar evolution
is under investigations.  



\acknowledgements
I acknowledge the suggestion by E. Schatzman and A. Baglin to 
undertake the development of this code; J. Provost and G. Berthomieu for 
advices, debugging
and constructive criticisms of the results; J. Christensen-Dalsgaard,
M. Gabriel, A. Noels, J. Reiter for stimulating and helpful discussions;  
C. van't Veer, Y. Lebreton, M.J. Goupil and my colleagues of GDR CNRS 131 for
various contributions to CESAM. I wish to express my thanks to the anonymous 
referee whose comments and remarks greatly helped to improve the presentation 
of this paper.

\begin{thebibliography}{}

\bibitem[1986]{al}
Alecian, G. 1986, A\&A 168, 204

\bibitem[1989]{ag}
Anders, E., Grevesse, N. 1989, Geochimica et Cosmochimica Acta, 53, 197

%\bibitem[1994]{ab} Antia, H.M., Basu, S. 1994, ApJ 426, 801

\bibitem[1969]{at}
Arnett, W.D. Truran, J.W. 1969, ApJ 157, 339

\bibitem[1994]{ap}
Audard, N., Provost, J. 1994, A\& A 282, 73

\bibitem[1993]{be}
Berthomieu, G., Provost, J., Morel, P., Lebreton, Y. 1993, A\&A 268, 775

\bibitem[1978]{db}
de\,Boor, C. 1978, A Practical Guide to Splines
Springer, third edition 1985

\bibitem[1973]{bs}
de\,Boor, C., Swartz, B. 1973, SIAM J. Numer. Anal. 10, 582


\bibitem[1993]{br}
Bressan, A., Fagotto, F., Bertelli, G., Chiosi, C. 1993, A\&AS 100, 647

\bibitem[1988]{cf}
Caughlan, G.R., Fowler, W.A. 1988,
Atomic Data and Nuclear Data Tables, 40, 284

\bibitem[1968]{cy}
Clayton, D.D. 1968, Principles of Stellar Evolution and
Nucleosynthesis, Mc Graw--Hill

\bibitem[1992]{cvz}
Charbonel, C., Vauclair, S., Zahn, J.P. 1992, A\&A 255, 191

\bibitem[1995]{ch}
Chmielewski, Y., Cayrel de Strobel, G., Cayrel, R., Lebreton, Y., Spite M. 1995,
A\&A 299, 809

\bibitem[1982]{cd0}
Christensen-Dalsgaard, J. 1982, MNRAS 199,735

\bibitem[1988]{cd1}
Christensen-Dalsgaard, J. 1988, Computational procedures for GONG solar model 
project, Astronomisk Institut, Aarhus Universitet

\bibitem[1991]{cd2}
Christensen-Dalsgaard, J. 1991,
Solar oscillations and the physics of the solar interior.
Challenges to Theories of the Structure of Moderate-Mass Stars, Eds.
D.O. Gough \& J. Toomre, Springer Verlag, p. 11

\bibitem[1968]{cy}
Clayton, D.D. 1968, Principles of Stellar Evolution and
Nucleosynthesis, Mc Graw--Hill


\bibitem[1987]{cdb}
Conte, S.D., de\,Boor, C. 1987, Elementary Numerical Analysis, McGraw-Hill Book 
Company, third edition

\bibitem[1989]{co}
Cox, A.N., Guzik, J.A., Kidman, R.B. 1989, ApJ 342, 1187

\bibitem[1971]{eg}
Eggleton, P. 1971, MNRAS 151, 351

\bibitem[1972]{eg1}
Eggleton, P. 1972, MNRAS 156, 361

\bibitem[1973]{eff}
Eggleton, P., Faulkner J., Flannery B.P. 1973, A\&A 23, 325

\bibitem[1989]{em}
Embarek, B. 1989, Stage de Maitrise d'In\-g\'e\-nie\-rie 
ma\-th\'e\-ma\-ti\-que, Universit\'e de Nice, Juin 1989

\bibitem[1990]{gab}
Gabriel, M. 1990, in Progress of seismology of the sun and
stars, eds. Osaki, Y., and Shibahashi, H., Springler, New York, p. 23

\bibitem[1988]{ga}
Gaches, J. 1988, Mod\'elisation de Courbes et surfaces Splines,
Ecole Nationale Sup\'erieure d'Ing\'enieurs de Constructions A\'eronautiques,
Toulouse

\bibitem[1971]{ge}
Gear, C.W. 1971, Numerical initial value problems in ordinary differential
equations, Prentice Hall


\bibitem[1993]{go}
Goupil, M.J., Michel, E., Lebreton, Y., Baglin, A. 1993, A\&A 268, 546

\bibitem[1993]{gr}
Grevesse, N., Noels, A. 1993, in Origin and Evolution of the Elements. Eds. N. 
Prantzos, E. Langioni-flam, M. Casse, Cambridge Univ. Press, p.14

\bibitem[1992]{gu}
Guenther, D.B., Demarque, P., Kim, Y.C., Pinsonneault, M.H. 1992, ApJ 387, 372

\bibitem[1995]{gm}
Guillot, T., Morel, P. 1995, A\&AS 109, 109

\bibitem[1991]{hw}
Hairer, E., Wanner, G. 1991, Solving Ordinary Differential 
Equations II, Springer-Verlag

\bibitem[1994]{ha}
Han, Z., Podsiadlowski, P., Eggleton, P.P. 1994, MNRAS 270, 121

\bibitem[1966]{hs}
H\" arm, R., Schwarzchild, M. 1966, ApJ 145, 1966

\bibitem[1959]{hl}
Henyey, L.G., Le Levier, R., Levee, R.D. 1959, ApJ 129, 2
 
\bibitem[1964]{hf}
Henyey, L.G., Forbes, J.E., Goult, N.L. 1964, ApJ 139, 306

\bibitem[1965]{hv}
Henyey, L.G., Vardya, M.S., Bodenheimer, P.L. 1965, ApJ 142, 841

\bibitem[1996]{ho}
Houdek, G., Rogl, J. 1996,  Bull. Astr. Soc. India, 24, 317

\bibitem[1975]{ib}
Iben, I. 1975, ApJ 196, 525

\bibitem[1968]{ks}
Kippenhahn, R., Weigert, A., Hofmeister, E. 1968,
Methods of Computational Physics, 7, 129

\bibitem[1991]{kw}
Kippenhahn, R., Weigert, A. 1991, Stellar Structure and Evolution,
Springer-Verlag

\bibitem[1991]{kz}
Kurucz, R.L. 1991,
Stellar Atmospheres: Beyond Classical Models
(ed. L.Crivellari, I. Hubeny, and D.G. Hummer), NATO ASI Series, Kluwer,
Dordrecht 1991

\bibitem [1995]{le}
Lebreton, Y., Michel, E., Goupil, M.J., Baglin, A., Fernandes, J. 1995,
Accurate parallaxes and stellar ages determinations
IAU  Symp.  166, La Haye, August 1994.  Dordrecht :  D.
Reidel; eds. E. Hoeg, P.K. Seidelmann, p.~135

\bibitem[1987]{mm}
Maeder, A., Meynet, G. 1987 A\&A 182, 243

\bibitem[1985]{mk}
Marchouk, G., Agochkov, V. 1985, M\'ethode des El\'ements Finis, Mir

\bibitem[1993]{mi}
Michaud, G., Proffitt, C.R. 1993, Inside the Stars, IAU colloquium 137, ASP
Conf. Ser., 40, Eds. W.W. Weiss and A. Baglin, p. 246

\bibitem[1978]{ml}
Mihalas, D. 1978, Stellar Atmosphere, 2d Edition, Freeman and Cie

\bibitem[1988]{mhd}
Mihalas, D., D\"appen, W., Hummer, D.G. 1988, ApJ, 331, 815

\bibitem[1989]{mp0}
Morel, P. 1989, Un code d'Evolution Stellaire Adaptatif et 
Modulaire,
Cours de Structure Interne, Ecole d'Aussois, Eds. A.M. Hubert \& E. Schatzman,
Publ. Obs. Paris section Meudon

\bibitem[1990]{mp1}
Morel, P., Provost, J., Berthomieu, G. 1990, Solar Physics 282, 73

\bibitem[1993]{mp2}
Morel, P., van't Veer, C., Berthomieu, G., Castelli, F., Cayrel, R., Lebreton, 
Y., Provost, J. 1993, Inside the Stars, IAU colloquium 137,
Eds. W.W. Weiss and A. Baglin, p. 57

\bibitem[1993a]{mp5}
Morel, P. 1993a, Utilisation et description du Code d'Evolution Stellaire 
CESAM, Pub. Obs. Nice

\bibitem[1993b]{mp3}
Morel, P. 1993b, Inside the Stars, IAU colloquium 137,
Eds. W.W. Weiss and A. Baglin, p. 445

\bibitem[1994]{mp4}
Morel, P., van't Veer, C., Provost, J. Berthomieu, G., Castelli, F., Cayrel, 
R., Lebreton, Y. 1994, A\&A 286, 91

\bibitem[1996]{mp6}
Morel, P., Schatzman, E. 1996, A\&A 310, 982

%\bibitem[1993]{neu} Neuforges, C. 1993, A\&A 268, 650

\bibitem[1995]{no}
Noels, A., Grevesse, N., Mazzitelli, I. 1995, Processes in 
Astrophysics, International meeting in honour of 
Evry Schatzman. Lecture notes in physics, Springler. Eds. I.W. Roxburg, J.L. 
Masnou, p.~203 

\bibitem[1969]{pa}
Paczynski, B. 1969, Acta Astron. 19, 1

\bibitem[1986]{pf}
Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, 
W.T. 1986, Numerical Receipes, Cambridge University Press, Cambridge

\bibitem[1991]{pm}
Proffitt, C.R., Michaud, G. 1991  Apj, 371, 584

\bibitem[1994]{qv}
Quarteroni, A., Valli, A. 1994, Numerical approximation of partial
differential equations, Springer-Verlag

\bibitem[1994]{re}
Reiter, J., Bulirsh, R., Pfleiderer, J. 1994, Astron. Nachr. 315, 205

\bibitem[1995]{re1}
Reiter, J., Walsh, L., Weiss, A. 1995, MNRAS 274, 899

\bibitem[1967]{rm}
Richtmyer, R.D., Morton, K.W. 1967, Finite Differences Methods for Initial 
Values Problems, 2d. Edt., Intersciences Publishers, New-York

\bibitem[1992]{ri}
Rogers, F.J., Iglesias, C. A. 1992, ApJS 79, 507

\bibitem[1995]{ro}
Rossignol, S. 1995, M\'emoire de DEA, Astronomie , Imagerie et Haute 
R\'esolution Angulaire, Universit\'e de Nice-Sophia Antipolis

\bibitem[1959]{s1}
Schatzman, E. 1959, Sur la perte de masse et les processus de freinage de la 
rotation, I.A.U. Symposium 10, ed. J. Greenstein, Ann. Astr. Supp. 8, 129

\bibitem[1981]{sm}
Schatzman, E., Maeder, A., Angrand, F., Glowinski, R. 1981, A\&A 96, 1

\bibitem[1990]{sp}
Schatzman, E., Praderie, F. 1990, Les Etoiles,
InterEditions/Editions du CNRS

\bibitem[1981]{sk}
Schumaker, L. 1981, Splines Functions: Basic Theory, John Whiley

\bibitem[1979]{sb}
Stoer, J., Bulirsch, R. 1979, Introduction to Numerical Analysis, 
Springer-Verlag.

\bibitem[1970]{st}
Strittmatter, P.A., Faulkner, J., Robertson, J.W., Faulkner, D.J. 1970,
ApJ 161, 369

\bibitem[1969]{wa}
Wagoner, R.V. 1969, ApJS 18, 247

\bibitem[1981]{wi}
Wilson, R.E. 1981, A\&A 99, 43

\bibitem[1991]{jpz}
Zahn, J.P. 1991, A\&A 252, 179

\end{thebibliography}

\input AA_article_app.tex

\end{document}

