Subroutine root_ml(f,a,b,eps,Roots,nint,n)
c====================================================================
c Program to find multiple roots of an equation f(x)=0
c using the Brute force method together with Bisectional method
c written by: Alex Godunov (October 2006)
c--------------------------------------------------------------------
c input ...
c f	- function which evaluates f(x) for any x in the interval a,b
c a	- left endpoint of initial interval
c b	- right endpoint of initial interval
c eps	- desired uncertainity of the root
c nint  - number of intervals to split (a-b)
c
c output ...
c Roots - roots of equation f(x)=0
c n     - number of roots (n<nint)
c
c The program split the interval (a-b) into (a-b)/nint intervals
c====================================================================
	Integer*4 nint, n
	Real*8 f,a,b,eps,Roots(nint)
	Real*8 xl, x0, xr, dx, ai, bi
	Integer*4 i,j, iflag, iter
	Parameter (iter=1000)
	External f
c* step for new interval
      dx = (b-a)/real(nint)
c* loop over new intervals
      n = 0
      do j=1,nint
         ai = a + real(j-1)*dx
         bi = ai + dx
         i = 0
         xl = ai
         xr = bi
c* bisectional method to the new interval if f(xl)*f(xr).le.0.0
         if(f(xl)*f(xr).le.0.0) then
            do while (abs(xr-xl).gt.eps)
               i=i+1
               x0=(xr+xl)/2
               if(f(xl)*f(x0).le.0.0) then
                  xr = x0
                  else
                  xl=x0
               end if
               if(i.ge.iter) exit
            end do
            n = n+1
            Roots(n)=x0
         end if
c* end of bisectional method
      end do
      return
      end