anchor.diff:               call variables(T, Smin, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e,&
                     tau_eff, P_rad, P_gaz,E_ff,Fz,f_min,Sigma_0, Omega_0,rs)
             endif
  +    
             dichotomy       = (Smin+Smax)/2.0d0
             call variables(T, dichotomy, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
                  K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f_dichotomy,Sigma_0, Omega_0,rs)
             j               = j + 1
  +    
          end do iteration
       endif
  +  
     end function dichotomy
     
     
  
positive.code:        !do j              = 1 , n_cell
         !  r = (param%rmax-rmin)/(n_cell-1)*(j-1) + rmin
         r              = 10*G*param%M/(c**2)     
         omega          = G*param%M/r**3
         write(number_of_cell,'(I5.5)') j
         fid = 20 + j
         fname = 'Temperature_Sigma_'//trim(number_of_cell)//'.dat'
         open(fid,file   = fname, status='unknown',action='readwrite') 
         !  write(fid,*)"T              Sigma"
          do i           = 1 , nb_it
            Smin        = 1.0d-8
            Smax        = 1.0d9
            temp        = (t_max-t_min)/(nb_it-1)*(i-1) + t_min
            sigma       = dichotomy(Smin, Smax, eps, temp, omega, sigma_0, Omega_0,rs)
            call variables(temp, sigma, omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
  
negative.code:     !-------------------------------------------------------------------------
      dichotomy             = (Smin+Smax)/2.
      j = 0
      
      call variables(T, Smin, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e,&
           tau_eff, P_rad, P_gaz,E_ff,Fz,f_min,Sigma_0, Omega_0,rs)
      call variables(T, Smax, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff, K_e, &
           tau_eff, P_rad, P_gaz,E_ff,Fz,f_max,Sigma_0, Omega_0,rs)
      call variables(T, dichotomy, Omega, H, rho, cs, nu, Q_plus, Q_minus, K_ff,&
           K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f_dichotomy,Sigma_0, Omega_0,rs)
      
      write(*,*)'fmin = ',f_min
      write(*,*)'fmax = ',f_max
      
      !call display_variables(T,Omega,r, sigma, H, rho, cs, nu, Q_plus, Q_minus,&
       !    K_ff, K_e, tau_eff, P_rad, P_gaz,E_ff,Fz,f)
      
      if ( f_max * f_min .gt. 0.) then
         write(*,*)'This function image does not switch its sign in this particular interval.'
      endif
  