PROGRAM BDF
!               ! Solves a system of ODEs using a backward
!               ! differentiation formula with Newton 
!               ! iteration
!--------------------------------------- J R Dormand, 1/96
! Modified: A Chisalita, 1/2006
! ------------------------------

  USE Stiff
  IMPLICIT NONE
  INTEGER :: neqs, i, j, p, pp1, pm1, nstep, its, ldiff, &
             interv
  INTEGER, PARAMETER :: itmax = 10
  INTEGER, ALLOCATABLE :: order(:)
  REAL (KIND = 8), ALLOCATABLE :: y(:, :), amat(:, :),&
   ident(:,:), yp(:), bdel(:, :), a(:),&
   alpha(:), beta(:), c(:), dy(:)
  REAL (KIND = 8) :: x, h, hrk, b, rb, xend, errfun, tol, &
                     error

character(260) fname, fname2, dummy*80
logical true_solution
real(8) x0, y0(neq)
integer iPoint



print '(a,\)', ' Input data file: '
read *, fname
open(1, file =fname, form ='formatted', status ='old')

read(1,'(a)') dummy
read(1,*) p	! BDF order
read(1,*) h, tol
read(1,*) x0, y0(1:neq)	! t, x0
read(1,*) xEnd
read(1,*) interv
read(1,*) true_solution

! Result file:
iPoint =index(fname,'.',back=.TRUE.)
fname2 =fname(1:iPoint-1)//'.rez'

open(2, file =fname2, form ='formatted', status ='unknown')
write(2,'(a)') '   t 		 X					Xprime				 Err(X)			Err(Xprime)'


  pm1 = p - 1; pp1 = p + 1
  ALLOCATE(a(p), alpha(0:p), beta(p))
!------Generate BDF coefficients for required order--------
  alpha = 0.0d0; beta = 0.0d0; alpha(0) = 1; rb=0.0d0
  DO i = 1, p
     DO j = i, 1, -1
        alpha(j) = alpha(j) - alpha(j-1)
        beta(j)  = beta(j) + alpha(j)/i
     END DO
     rb = rb + 1.0d0/i; b = 1.0d0/rb
     a = -b*beta
  END DO
  neqs = neq
!  tol = 1.0d-6             ! Tolerance for Newton iteration
  ALLOCATE( y(neqs, pp1), amat(neqs, neqs),  &
    yp(neqs), order(neqs),c(neqs), bdel(neqs, pp1),&
    dy(neqs), ident(neqs,neqs)) 
  ident = 0.0d0;  DO i = 1, neq; ident(i,i) = 1.0d0; END DO 

x =x0; y(:,2) =y0;
  nstep = (xend - x)/h + 0.5

!--------------------- Starter values ---------------------
  DO i = 3, pp1
     hrk = h
     y(:, i) = RKstart(y(:, i-1), hrk)  ! Runge-Kutta 
  END DO
  ldiff = 2   ! Extrapolation is lower order for 1st step
!--------------------- Start integration ------------------
  DO j = p, nstep                         ! Loop over steps
! --- Predictor is polynomial extrapolation formula 
! ---------------------------------- to avoid f evaluations
    bdel = y
    yp = y(:, pp1)
    DO i = 1, pp1 - ldiff  
! --------------------------------------Compute differences 
       bdel(:, 1:pp1 - i) = bdel(:, 2: pp1-i+1)&
                                    - bdel(:, 1: pp1-i)
! -----------------Accumulate interpolant to form predictor
       yp = yp + bdel(:, pp1 - i)
    END DO
    ldiff = 1    ! Predictor order same as BDF after step 1
    c = 0
    DO i = 0, pm1;  c = c + a(i+1)*y(:, pp1 - i);  END DO
!----------------------------------------------------------
! -----Correct with BDF formula using Newton iteration ----                   
    its = 0
    x = x + h
    error = 1.0d0
    DO While(error > tol)       ! Commence Newton iteration
      amat = ident - h*b*Jac(x, yp)          ! Use Jacobian
      dy = -yp + h*b*Fcn(x, yp) + c          ! f evaluation
      errfun = MAXVAL(Abs(dy))  
      CALL LU(amat, neq, order)         ! Solve Newton equns 
      CALL FBSubst(amat, dy, neq, order)   ! by LU factoring
      error = MAXVAL(Abs(dy)) + errfun     ! error statistic
      yp = yp + dy;   its = its + 1;       ! Update solution
      IF(its > itmax) THEN
         PRINT*, 'Maximum iterations exceeded'
         STOP
      END IF
    END DO
    DO i = 1, p; y(:, i) = y(:, i+1); END DO ! Update
    y(:, pp1) = yp                          ! starter values
    IF(MOD(j, interv) == 0) THEN
       CALL Output(2, x, y(:, pp1), true_solution)
    ENDIF

	if(j==nstep) then
	write(1,'(/,a)') ' Results at tEnd:'
	write(1,'(a)') '  tEnd		 X					Xprime				Err(X)			Err(Xprime) '
       CALL Output(1, x, y(:, pp1), true_solution)
	endif
  END DO
  CONTAINS
!----------------------------------------------------------
  FUNCTION RKstart(y, h)   ! Carries out RK integration  to
                           ! hit x + h. Several steps are
                           ! computed if necessary according
                           ! to local tolerance 1.0e-6.
!-----------------------------------------J R Dormand, 1/96
!
  IMPLICIT NONE
  INTEGER ::  ir, istep
  LOGICAL :: hit
  REAL (KIND = 8) :: y(:), h, Rkstart(SIZE(y))
  REAL (KIND = 8) :: letol, errest, xhit, beta
  REAL (KIND = KIND(1.0d0)), DIMENSION(2: 5, 4),&
       PARAMETER :: a =  RESHAPE((/0.2d0, 0.0d0, 1.2d0,&
          -2.125d0, 0.0d0, 0.4d0, -2.4d0, 5.0d0, 0.0d0,&
    0.0d0, 2.0d0, -2.5d0, 0.0d0, 0.0d0, 0.0d0, 0.625d0/),&
    (/4, 4/))
  REAL (KIND= KIND(1.0d0)), DIMENSION(5), PARAMETER :: b =&
   (/ 13d0/96d0, 0.0d0, 25d0/48d0, 25d0/96d0, 1d0/12d0 /),&
      c = (/0.0d0, 0.2d0, 0.4d0, 0.8d0, 1.0d0/), &
   d = (/ 3d0/192d0, 0.0d0, -5d0/96d0, 15d0/192d0,&
        -1d0/24d0 /)
  REAL (KIND =8), ALLOCATABLE :: ff(:,:), w(:), ys(:)
  ALLOCATE(ff(5, neqs), w(neqs), ys(neqs))
    ys = y
    xhit = x + h; hit = .TRUE.
    letol = 1.0d-6
    ir = 0; istep = 0
    ff(1, :) = Fcn(x, ys)
    DO
      DO j = 2, 5
         w = ys + h*MATMUL(a(j, 1: j - 1), ff(1: j - 1, :))
         ff(j, :) = Fcn(x + c(j)*h, w)
      END DO
      errest = h*MAXVAL(ABS(MATMUL(d, ff)))
      IF ( errest < letol) THEN
         ys = ys + h*MATMUL(b, ff)
         istep = istep + 1
         IF(hit) THEN
           x = xhit
           PRINT '(1x, a, i3)', 'Steps =', istep
           RKstart = ys
           RETURN
         ENDIF
         x = x + h
         beta = errest/letol
         h = MIN(0.9d0*h/beta**0.25, 2*h)
         IF(xhit < x + h) THEN
           h = xhit - x
           hit = .TRUE.
         END IF
         ff(1, :) = Fcn(x, ys)
      ELSE
         hit = .FALSE. 
         beta = errest/letol
         h = MAX(0.9d0*h/beta**0.25, h/10)
         ir = ir + 1
         PRINT '(1x,a,i3)', 'Rejects in RKstart =', ir
      ENDIF
    END DO
  END FUNCTION Rkstart   
 !
 END PROGRAM BDF
