recsAuxiliary calculates auxiliary variables not included in the core model
 
  INTERP = recsAuxiliary(MODEL,INTERP,S,X,Z)
  For the definition of MODEL and INTERP structure, see recsSolveREE
  documentation. recsAuxiliary returns the interpolation structure INTERP, in
  which the coefficients matrices cxa (for auxiliary variables) and cza (for
  auxiliary expectations) have been added as new fields.
 
  INTERP = recsAuxiliary(MODEL,INTERP,S,X,Z,XA) uses XA as first guess for the
  value of the auxiliary variables on the grid. This is only useful when
  auxiliary variables are function of auxiliary expectations. In this case, it
  can reduce a lot the time to reach the solution to start from a good first
  guess.
 
  INTERP = recsAuxiliary(MODEL,INTERP,S,X,Z,XA,OPTIONS) solves for auxiliary
  variables with the parameters defined by the structure OPTIONS. The fields of
  the structure are
     eqsolver         : 'fsolve', 'krylov' (default), 'lmmcp', 'ncpsolve',
                        'path' or 'sa'
     extrapolate      : 1 or 2 if extrapolation is allowed outside the
                        interpolation space, 0 or -1 to forbid it (default: 1).
     eqsolveroptions  : options structure to be passed to eqsolver
 
  [INTERP,XA] = recsAuxiliary(MODEL,INTERP,S,X,Z,...) returns the value of the
  auxiliary variables on the grid.
 
  [INTERP,XA,ZA] = recsAuxiliary(MODEL,INTERP,S,X,Z,...) returns the value of
  the auxiliary expectations on the grid.