anchor.diff:         iteration:do while ( dabs( Smax - Smin ) .ge. eps .and. j .lt. 10000)
   
   
  -    call variables(T, Smin, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e,&
  -         tau_eff, E_ff,Fz,f_min, optical_depth)
  -
  -    call variables(T, Smax, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e, &
  -         tau_eff, E_ff,Fz,f_max, optical_depth)
  -
  -    call variables(T, S_center, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
  -         K_e, tau_eff, E_ff,Fz,f_center, optical_depth)
  +    call variables(T, Smin, Omega, f_min, optical_depth)
  +
  +    call variables(T, Smax, Omega, f_max, optical_depth)
  +
  +    call variables(T, S_center, Omega, f_center, optical_depth)
   
   
             if(f_min * f_center .gt. 0.) then
  
positive.code:   !Subroutine in order to compute variables H, rho, cs, nu, Q_plus, Q_minus,
    !K_ff, K_e, tau_eff, E_ff,Fz given T, Sigma and Omega
    !------------------------------------------------------------------------
    subroutine variables(T, Sigma, Omega, H, rho, cs, nu, Q_plus, Q_minus,&
         K_ff, K_e, tau_eff, E_ff,Fz,f, optical_depth)
      implicit none
  
      real(kind = x_precision),intent(in)                      :: T,Sigma,Omega
      integer, intent(in)                                      :: optical_depth
      real(kind = x_precision)                                 :: coeff_a=0.,coeff_b=0.,coeff_c=0.
  
      real(kind = x_precision),intent(out)                     :: H
      real(kind = x_precision),intent(out)                     :: rho
      real(kind = x_precision),intent(out)                     :: cs
      real(kind = x_precision),intent(out)                     :: nu
      real(kind = x_precision),intent(out)                     :: Q_plus
      real(kind = x_precision),intent(out)                     :: K_ff
      real(kind = x_precision),intent(out)                     :: K_e
      real(kind = x_precision),intent(out)                     :: E_ff
      real(kind = x_precision),intent(out)                     :: tau_eff
      real(kind = x_precision),intent(out)                     :: Fz
      real(kind = x_precision),intent(out)                     :: Q_minus
      real(kind = x_precision),intent(out)                     :: f
      !------------------------------------------------------------------------
  
      coeff_a              = (Omega * state_0%Omega_0)**2 * Sigma * state_0%S_0 / 2._x_precision
  
negative.code: 
            sigma       = dichotomy(Smin, Smax, eps, temp, omega, optical_depth)
  
            call variables(temp, sigma, omega, f, optical_depth)
  
            temp_t_1(i)   =  temp
            sigma_t_1(i)  = sigma
  