z1 = dble(fp) + dble(fpa) - 3.d0*dble(fa-f)/dble(ta-t)

         discri = z1*z1-dble(fp)*dble(fpa)

      subroutine CUBIC( t, f, fp, ta, fa, fpa, tlower, tupper )

         discri = discri*dble(fpa)

         den  = sngl(z1+b+dble(fpa))

      double precision    t, f, fp, ta, fa, fpa, tlower, tupper

      fpa    = dotdg

         call CUBIC( tact, fnew, fp, ta, fa, fpa, left, right )

         call CUBIC( tact, fnew, fp, ta, fa, fpa, left, right )

         fpa = fp

      double precision    fa, fpa, fp