!Prediction-Correction algorithm routine for a single first order 
!differential equation.
!s = steplength, c = correction term.

subroutine PC_general(f, x1, y1, s, c, x2, y2, err, d)
implicit none

     integer :: i
     integer, intent(in) :: d
     real :: y2_ini(d), y2_prime_ini(d), y1_prime(d), diff(d)
     real, intent(in) :: x1, y1(d), s 
     real, intent(out) :: x2, y2(d), err(d)
     real, intent(inout)::  c(d)
     real, external :: f


!Initial step, (x0,y0) -> (x1,y1) (Runge-Kutta), must be 
!taken by driver routine.

x2 = x1 + s
do i = 1, d
   y1_prime(i) = f(x1, y1(1), y1(2), i)
enddo

!Prediction step:
do i = 1, d
   y2_ini(i) = y1(i) + y1_prime(i)*s + c(i)*s**2
enddo

do i = 1, d
   c(i) = 0.5*(f(x2, y2_ini(1), y2_ini(2), i) - y1_prime(i))
enddo

!Correction step with corrected c:
do i = 1, d
   y2(i) = y1(i) + y1_prime(i)*s + c(i)*s
enddo

do i = 1, d
   diff(i) = y2_ini(i) - y2(i)
   err(i) = ABS(diff(i))
enddo

end subroutine PC_general