function aitken (f, p0)
  implicit none

  double precision f, p0, aitken
  external f

  double precision p1, p2

  p1 = f(p0)
  p2 = f(p1)
  aitken = p0 - (p1 - p0)**2 / (p2 - (2.0d0 * p1) + p0)
end

subroutine steff (f, pinit, tol, mxiter, ierr, pfinal)
  implicit none

  double precision f, pinit, tol, pfinal
  integer mxiter, ierr
  external f

  double precision p0, p, aitken
  integer iter
  external aitken

  p0 = pinit
  p = aitken (f, p0)
  iter = 1
  while (abs (p - p0) > tol && iter < mxiter)
    {
      p0 = p
      p = aitken (f, p0)
      iter = iter + 1
    }
  if (abs (p - p0) > tol)
    ierr = -1
  else
    {
      ierr = 0
      pfinal = p
    }
end

function dcstlj (c0, c1, c2, t)
  implicit none

  double precision c0, c1, c2, t, dcstlj

  double precision s, c01, c12, c012

  s = 1.0d0 - t
  c01 = (s * c0) + (t * c1)
  c12 = (s * c1) + (t * c2)
  c012 = (s * c01) + (t * c12)
  dcstlj = c012
end

function x (t)
  implicit none
  double precision x, t, dcstlj
  external dcstlj
  x = dcstlj (2.0d0, -8.0d0, 2.0d0, t)
end

function y (t)
  implicit none
  double precision y, t, dcstlj
  external dcstlj
  y = dcstlj (1.0d0, 2.0d0, 3.0d0, t)
end

function impleq (x, y)
  implicit none
  double precision x, y, impleq
  impleq = (5.0d0 * x * x) + y - 5.0d0
end

function f (t)
  implicit none
  double precision f, t, x, y, impleq
  external x, y, impleq
  f = impleq (x(t), y(t)) + t
end

program RCstef
  implicit none

  double precision x, y, impleq, f
  external x, y, impleq, f, steff

  double precision t0, t
  integer i, ierr

  t0 = 0.0d0
  for (i = 0; i != 11; i = i + 1)
    {
      call steff (f, t0, 0.00000001d0, 1000, ierr, t)
      if (ierr < 0)
        write (*,*) "t0 = ", t0, " : no answer"
      else if (abs (impleq (x(t), y(t))) <= 0.000001)
        write (*,*) "t0 = ", t0, " : intersection at (", _
                    x(t), ", ", y(t), ")"
      else
        write (*,*) "t0 = ", t0, " : spurious solution"
      t0 = t0 + 0.1d0
    }
end