C_______________________________________________________________________________
C
C     Este programa calcula la raiz de una funcion fprueba(x)
C     mediante el metodo de biseccion
C_______________________________________________________________________________

      implicit real*8 (a-h,o-z)
      external fprueba        !nombre de la funcion cuya raiz se pretende hallar 

	write(6,'(/a,$)') '  Extremo izquierdo del intervalo -> x      = '
      read(5,*) x
	write(6,'( a,$)') '  Extremo derecho   del intervalo -> z      = '
      read(5,*) z
	write(6,'( a,$)') '  Err. Relativo                   -> r      = '
      read(5,*) r
	write(6,'( a,$)') '  Err. Absoluto                   -> e      = '
      read(5,*) e
	write(6,'( a,$)') '  No. Max. Iter.                  -> mxiter = '
      read(5,*) mxiter
      write(6,'()')

      call Biseccion(fprueba,x,z,a,da,r,e,mxiter,niter)

      if (niter.lt.0) then
        write(6,'(a)') '  Error: el intervalo es inadecuado.'
      else if (niter.le.mxiter) then
        write(6,100) niter,x,fprueba(x),z,fprueba(z),a,da,fprueba(a)
  100   format(///' SOLUCION: '//
     .         '  niter = ',i6,'  -> x  = ',d24.16,'     ',24x,
     .                         ';  f(x) = ',d24.16/
     .         '          ',6x,'  -> z  = ',d24.16,'     ',24x,
     .                         ';  f(z) = ',d24.16/
     .         '          ',6x,'  -> a  = ',d24.16,' +/- ',d24.16,
     .                         ';  f(a) = ',d24.16)
      else
        write(6,'(a)') '  Error: el metodo no ha convergido.'
      endif

      read(5,'()')

      end


c______________________________________________fprueba(x)=cos(x)-x____________
c     Función que utilizamos como ejemplo

      real*8 function fprueba(x)
      implicit real*8 (a-h,o-z)

      fprueba=cos(x)-x

      return
      end


C_______________________________________________________________________________
C
C     Subrutina Biseccion(f,x,z,a,da,r,e,mxiter,niter)
C
C        f      -> entrada: funcion cuya raíz se pretende hallar
C        x      -> entrada: extremo izquierdo del intervalo inicial
C               -> salida : extremo izquierdo del intervalo final
C        z      -> entrada: extremo derecho   del intervalo inicial
C               -> salida : extremo derecho   del intervalo final
C        a      -> salida : solucion aproximada de la ecuación f(a)=0
C        da     -> salida : maximo error absoluto de la solucion aproximada 
C        r      -> entrada: tolerancia de la sol. aprox. en error relativo 
C        e      -> salida : tolerancia de la sol. aprox. en error absoluto 
C        mxiter -> entrada: no. maximo de iteraciones permitidas
C        niter  -> salida : no. de iteraciones realizadas
C                           [niter < 0      => intervalo inicial inadecuado]
C                           [niter > mxiter => el proceso no ha convergido]
C
C     Se utiliza el metodo de biseccion
C_______________________________________________________________________________

      Subroutine Biseccion(f,x,z,a,da,r,e,mxiter,niter)
      implicit real*8 (a-h,o-z)
      data cero /0.d+00/, dos /2.d+00/

      fx=f(x)
      fz=f(z)
      p=fx*fz
      if (p.gt.cero) then
        niter=-1   ! se indica asi que el intervalo inicial es inadecuado
        return
      else
        niter=0
      endif

      do while (niter.lt.mxiter)
        a=(x+z)/dos
        fa=f(a)
        da=(z-x)/dos
        ea=abs(a)*r
        if (e.gt.ea) ea=e
        write(6,100) niter,x,fx,z,fz,a,da,fa
  100   format(/
     .         '  niter = ',i6,'  -> x  = ',d24.16,'     ',24x,
     .                         ';  f(x) = ',d24.16/
     .         '          ',6x,'  -> z  = ',d24.16,'     ',24x,
     .                         ';  f(z) = ',d24.16/
     .         '          ',6x,'  -> a  = ',d24.16,' +/- ',d24.16,
     .                         ';  f(a) = ',d24.16)
        p=fx*fa
        if (da.lt.ea.or.p.eq.cero) return
        niter=niter+1
        if (p.lt.cero) then
          z=a
          fz=fa
        else
          x=a
          fx=fa
        endif
      enddo

      niter=mxiter+1   ! se indica asi que el proceso no ha convergido
      return

      end



