C==============================================================
C A simple algorithm for computing zero of a nonlinear function
C of a varioable in a given interval [a,b]
C
C Neculai Andrei
C March 28, 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

	f = x*dexp(x)+2.d0*dexp(x)-1.d0

	return
	end


c=========================================== Main program
c                                           ==============
	real*8 a,b,c, x,f
	real*8 fa,fb,fc,fx
	real*8 epsm
	integer iprint

	epsm=0.0000000001d0
	iprint = 0

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

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

	iter=0

	call func(a,fa)
	call func(b,fb)
	write(4,50) a,fa
50	format(4x,'a=',e20.13,5x,'fa=',e20.13)
	write(4,51) b,fb
51	format(4x,'b=',e20.13,5x,'fb=',e20.13,/)

	  if(fa*fb .gt. 0.d0) then
	write(4,100) a,b
100	format(4x,'Aparently in interval [',e20.13,',',e20.13,'] NO Zero')
	write(4,101)
101   format(4x,'Modify the interval by changing the bounds a and/or b')
	go to 998
	  end if

c------------------------------------------------------------
10	continue
	
	if(dabs(fa) .le. epsm .or. dabs(fb) .le. epsm) go to 999

	iter=iter+1

	if(fa .le. 0.d0) then
	  
	  c = a + (b-a)/2.d0
	  call func(c,fc)

	  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

	  c = a + (b-a)/2.d0
	  call func(c,fc)

	  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,113) b-a
113		  format(4x,'fa>0. b-a=',e20.13)
          end if
	    go to 10
	  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 10
	  end if
	end if

999	continue

	if(dabs(fa) .le. epsm) x = a
	if(dabs(fb) .le. epsm) x = b
	call func(x,fx)

	write(4,120) x,fx
120	format(2x,'Solution: ',e20.13,4x,'Function value: ',e20.13)
	write(4,121) iter
121	format(2x,'Number of iterations; ',i7)

998	continue

	stop
	end
c========================================================================
c Last line