C==============================================================
C A simple algorithm for computing all the zeros of a nonlinear 
C function of a variable in a given interval [a,b]. b>a.
C
C The algorithms tries to determine all the zeros of a nonlinear 
C function in a given interval by discretizing of this interval
C and by find the intervals where the function changes its sign.
C
C Every subinterval is reduced by its halving for finding the
C zero of the function under accuracy epsm.
C
C Neculai Andrei
C April 16, 1975
C==============================================================


	subroutine func(x,f)
	real*8 x,f

c	f = x**2 - 4.d0
c	f = x**2 - 8.d0*x - 9.d0
c	f = 4.d0*x**3 - 3.d0*x - 1.d0
c	f = x**3 + 4.d0*x**2 - 4.d0*x - 16.d0
c	f = 1.d0/((x-0.3d0)**2+0.01) + 1.d0/((x-0.9d0)**2+0.04) -6.d0
c	f = x**4 - 5.d0*x**2 + 4.d0
c	f = x**6 - 14.d0*x**4 + 49.d0*x**2 - 36.d0
c	f = x**7 + 7.d0*x**6 - 14.d0*x**5 - 98.d0*x**4 +
c     +    49.d0*x**3 +343.d0*x**2 - 36.d0*x - 252.d0
c	f = x**3 - 2.d0*x**2 + 1.d0
c	f = cos(x) - x**2
c	f = x*dexp(x) - 7.d0	
c	f = x - cos(x)
c	f = cos(x)
c	f = sin(x)
c	f = cos(x**2)
	f = x**3 + 3.d0*x**2 - 18.d0*x

	return
	end


c=========================================== Main program
c                                           ==============
	real*8  a, b, c, x
	real*8 fa,fb,fc,fx
	real*8 h
	real*8 epsm
	integer iprint
	integer i, ni, iint
	integer iter, itertot
	integer fgcnt, fgcntot

	epsm=0.000000001d0
	iprint = 0

C Searching interval
c ==================
	a = -7.d0
	b =  7.d0

      open(unit=4,file='zeros.out',status='unknown')

	itertot = 0
	fgcntot = 0

	write(4,1)
1	format(4x,'All zeros of a nonlinear function')
	write(4,2) a,b
2	format(4x,'in interval [',e20.13,','e20.13,' ]')
	write(4,3)
3	format(4x,56('='),/)

C-----------------------------------------------------
	call func(a, fa)
	fgcnt = fgcnt + 1


C Take a discretization of the interval [a,b]
C ===========================================

	ni = 200
	h = (b-a)/float(ni)


c Start iteratuions
c =================
	iint= 0
	i   = 1

400	continue
	iter  = 0
	fgcnt = 0

	if(i .gt. ni) go to 998

30	  b = a + h
	  iter = iter+1
	  call func(b, fb)
	  fgcnt = fgcnt + 1

	  if(fa*fb .gt. 0.d0) then
	    a  = b
	    fa = fb
	    i = i+1
	    if(i .gt. ni) go to 998
	    go to 30
	  else
c			write(4,104)
c104			format(4x,'fa*fb le 0.d0: Change of sign------')

		iint = iint+1
	    if(fa .le. 0.d0) then 
10		  continue
	    if(dabs(fa) .le. epsm .or. dabs(fb) .le. epsm) go to 999
	    iter=iter+1
	      c = a + (b-a)/2.d0
	      call func(c,fc)
	      fgcnt = fgcnt + 1
	      if(fc .ge. 0.d0) then
	        b = c
	        fb = fc
	        if(iprint .eq. 1) then
	          write(4,99) iter
99		      format(4x,'iter='i5)
		      write(4,110) a,b,c
110		      format(4x,' a='e20.13,4x,' b=',e20.13,4x,' c=',e20.13)
		      write(4,111) fa,fb,fc
111		      format(4x,'fa='e20.13,4x,'fb=',e20.13,4x,'fc=',e20.13)
		      write(4,112) b-a
112	          format(4x,'fa<0. b-a=',e20.13)
              end if
	        go to 10
	      else
	        a = c
	        fa = fc
	        if(iprint .eq. 1) then
	          write(4,99) iter
		      write(4,110) a,b,c
		      write(4,111) fa,fb,fc
		      write(4,112) b-a
              end if
	        go to 10
	      end if	  			
	    end if

	    if(fa. ge. 0.d0) then
20		  continue
	      if(dabs(fa) .le. epsm .or. dabs(fb) .le. epsm) go to 999
	      iter=iter+1
	      c = a + (b-a)/2.d0
	      call func(c,fc)
	      fgcnt = fgcnt + 1
	      if(fc .ge. 0.d0) then
	        a = c
	        fa = fc
	        if(iprint .eq. 1) then
	          write(4,99) iter
		      write(4,110) a,b,c
		      write(4,111) fa,fb,fc
		      write(4,112) b-a
              end if
	        go to 20
	      else
	        b  = c
	        fb = fc
	        if(iprint .eq. 1) then
	          write(4,99) iter
		      write(4,110) a,b,c
		      write(4,111) fa,fb,fc
		      write(4,112) b-a
              end if
	        go to 20
	      end if
	    end if
	  end if
c	end do
	    

c------------------------------------------------------------

	
999	continue

	if(dabs(fa) .le. epsm .or. dabs(fb) .le. epsm) x=a
	call func(x, fx)
	fgcnt = fgcnt + 1

	write(4,119) iint
119	format(2x,'Zero #',i2)
	write(4,120) x, fx
120	format(2x,'ZERO: ',e20.13,4x,'Function value: ',e20.13)
	write(4,121) iter
121	format(2x,'Number of iterations:          ',i7)
	write(4,122) fgcnt
122	format(2x,'Number of function evaluations:',i7,/)

	a  = x + h
	call func(a,fa)
c	fgcnt = fgcnt + 1

	i = i+1

	itertot = itertot + iter
	fgcntot = fgcntot + fgcnt

	go to 400

998	continue

	write(4,130)
130	format(2x,39('='))
	write(4,131) itertot
131	FORMAT(2X,'TOTAL # OF iterations:          ',i6)
	write(4,132) fgcntot
132	format(2x,'TOTAL # of function evaluations:',i6)
	stop
	end
c========================================================================
c Last line