dp86co Subroutine

    subroutine dp86co(me,x,y,xend,hmax,h,rtol,atol,itol,iprint, &
                        iout,idid,nmax,nstiff,safe, &
                        beta,fac1,fac2, &
                        nfcn,nstep,naccpt,nrejct)

    implicit none

    class(dop853_class),intent(inout)   :: me
    real(wp),intent(inout)              :: x
    real(wp),dimension(:),intent(inout) :: y
    real(wp),intent(in)                 :: xend
    real(wp),intent(inout)              :: hmax
    real(wp),intent(inout)              :: h
    real(wp),dimension(:),intent(in)    :: rtol
    real(wp),dimension(:),intent(in)    :: atol
    integer,intent(in)                  :: itol
    integer,intent(in)                  :: iprint
    integer,intent(in)                  :: iout
    integer,intent(out)                 :: idid
    integer,intent(in)                  :: nmax
    integer,intent(in)                  :: nstiff
    real(wp),intent(in)                 :: safe
    real(wp),intent(in)                 :: beta
    real(wp),intent(in)                 :: fac1
    real(wp),intent(in)                 :: fac2
    integer,intent(inout)               :: nfcn
    integer,intent(inout)               :: nstep
    integer,intent(inout)               :: naccpt
    integer,intent(inout)               :: nrejct

    real(wp),dimension(me%n) :: y1,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10
    real(wp) :: atoli,bspl,deno,err,err2,erri,expo1,fac,fac11,&
                facc1,facc2,facold,hlamb,hnew,posneg,rtoli,&
                sk,stden,stnum,xout,xph,ydiff
    integer :: i,iasti,iord,irtrn,j,nonsti,nrd
    logical :: reject,last,event,abort

    ! initializations
    nrd = me%nrdens
    facold = 1.0e-4_wp
    expo1 = 1.0_wp/8.0_wp - beta*0.2_wp
    facc1 = 1.0_wp/fac1
    facc2 = 1.0_wp/fac2
    posneg = sign(1.0_wp,xend-x)
    irtrn = 0     ! these were not initialized
    nonsti = 0    ! in the original code
    xout = 0.0_wp !

    ! initial preparations
    atoli = atol(1)
    rtoli = rtol(1)
    last = .false.
    hlamb = 0.0_wp
    iasti = 0
    call me%fcn(x,y,k1)
    hmax = abs(hmax)
    iord = 8
    if ( h==0.0_wp ) then
        h = me%hinit(x,y,posneg,k1,iord,hmax,atol,rtol,itol)
    end if
    nfcn = nfcn + 2
    reject = .false.
    me%xold = x
    if ( iout/=0 ) then
        irtrn = 1
        me%hout = 1.0_wp
        call me%solout(naccpt+1,me%xold,x,y,irtrn,xout)
        abort = ( irtrn<0 )
    else
        abort = .false.
    end if

    if (.not. abort) then
        do
            ! basic integration step
            if ( nstep>nmax ) then
                if ( iprint/=0 ) &
                        write (iprint,'(A,E18.4)') ' exit of dop853 at x=', x
                if ( iprint/=0 ) &
                        write (iprint,*) ' more than nmax =' , nmax , 'steps are needed'
                idid = -2
                return
            elseif ( 0.1_wp*abs(h)<=abs(x)*uround ) then
                if ( iprint/=0 ) &
                        write (iprint,'(A,E18.4)') ' exit of dop853 at x=', x
                if ( iprint/=0 ) &
                        write (iprint,*) ' step size too small, h=' , h
                idid = -3
                return
            else
                if ( (x+1.01_wp*h-xend)*posneg>0.0_wp ) then
                    h = xend - x
                    last = .true.
                end if
                nstep = nstep + 1
                ! the twelve stages
                if ( irtrn>=2 ) call me%fcn(x,y,k1)
                y1 = y + h*a21*k1
                call me%fcn(x+c2*h,y1,k2)
                y1 = y + h*(a31*k1+a32*k2)
                call me%fcn(x+c3*h,y1,k3)
                y1 = y + h*(a41*k1+a43*k3)
                call me%fcn(x+c4*h,y1,k4)
                y1 = y + h*(a51*k1+a53*k3+a54*k4)
                call me%fcn(x+c5*h,y1,k5)
                y1 = y + h*(a61*k1+a64*k4+a65*k5)
                call me%fcn(x+c6*h,y1,k6)
                y1 = y + h*(a71*k1+a74*k4+a75*k5+a76*k6)
                call me%fcn(x+c7*h,y1,k7)
                y1 = y + h*(a81*k1+a84*k4+a85*k5+a86*k6+a87*k7)
                call me%fcn(x+c8*h,y1,k8)
                y1 = y + h*(a91*k1+a94*k4+a95*k5+a96*k6+a97*k7+a98*k8)
                call me%fcn(x+c9*h,y1,k9)
                y1 = y + h*(a101*k1+a104*k4+a105*k5+a106*k6+a107*k7+&
                            a108*k8+a109*k9)
                call me%fcn(x+c10*h,y1,k10)
                y1 = y + h*(a111*k1+a114*k4+a115*k5+a116*k6+a117*k7+&
                            a118*k8+a119*k9+a1110*k10)
                call me%fcn(x+c11*h,y1,k2)
                xph = x + h
                y1 = y + h*(a121*k1+a124*k4+a125*k5+a126*k6+a127*k7+&
                            a128*k8+a129*k9+a1210*k10+a1211*k2)
                call me%fcn(xph,y1,k3)
                nfcn = nfcn + 11
                k4 = b1*k1+b6*k6+b7*k7+b8*k8+b9*k9+b10*k10+b11*k2+b12*k3
                k5 = y + h*k4
                ! error estimation
                err = 0.0_wp
                err2 = 0.0_wp
                if ( itol==0 ) then
                    do i = 1 , me%n
                        sk = atoli + rtoli*max(abs(y(i)),abs(k5(i)))
                        erri = k4(i) - bhh1*k1(i) - bhh2*k9(i) - bhh3*k3(i)
                        err2 = err2 + (erri/sk)**2
                        erri = er1*k1(i) + er6*k6(i) + er7*k7(i) + er8*k8(i) &
                                + er9*k9(i) + er10*k10(i) + er11*k2(i) &
                                + er12*k3(i)
                        err = err + (erri/sk)**2
                    end do
                else
                    do i = 1 , me%n
                        sk = atol(i) + rtol(i)*max(abs(y(i)),abs(k5(i)))
                        erri = k4(i) - bhh1*k1(i) - bhh2*k9(i) - bhh3*k3(i)
                        err2 = err2 + (erri/sk)**2
                        erri = er1*k1(i) + er6*k6(i) + er7*k7(i) + er8*k8(i) &
                                + er9*k9(i) + er10*k10(i) + er11*k2(i) &
                                + er12*k3(i)
                        err = err + (erri/sk)**2
                    end do
                end if
                deno = err + 0.01_wp*err2
                if ( deno<=0.0_wp ) deno = 1.0_wp
                err = abs(h)*err*sqrt(1.0_wp/(me%n*deno))
                ! computation of hnew
                fac11 = err**expo1
                ! lund-stabilization
                fac = fac11/facold**beta
                ! we require  fac1 <= hnew/h <= fac2
                fac = max(facc2,min(facc1,fac/safe))
                hnew = h/fac
                if ( err<=1.0_wp ) then
                    ! step is accepted
                    facold = max(err,1.0e-4_wp)
                    naccpt = naccpt + 1
                    call me%fcn(xph,k5,k4)
                    nfcn = nfcn + 1
                    ! stiffness detection
                    if ( mod(naccpt,nstiff)==0 .or. iasti>0 ) then
                        stnum = 0.0_wp
                        stden = 0.0_wp
                        do i = 1 , me%n
                            stnum = stnum + (k4(i)-k3(i))**2
                            stden = stden + (k5(i)-y1(i))**2
                        end do
                        if ( stden>0.0_wp ) hlamb = abs(h)*sqrt(stnum/stden)
                        if ( hlamb>6.1_wp ) then
                            nonsti = 0
                            iasti = iasti + 1
                            if ( iasti==15 ) then
                                if ( iprint/=0 ) &
                                        write (iprint,*) &
                                        ' the problem seems to become stiff at x = ', x
                                if ( iprint==0 ) then
                                    idid = -4      ! fail exit
                                    return
                                end if
                            end if
                        else
                            nonsti = nonsti + 1
                            if ( nonsti==6 ) iasti = 0
                        end if
                    end if
                    ! final preparation for dense output
                    event = (iout==3) .and. (xout<=xph)
                    if ( iout==2 .or. event ) then
                        ! save the first function evaluations
                        do j = 1 , nrd
                            i = me%icomp(j)
                            me%cont(j) = y(i)
                            ydiff = k5(i) - y(i)
                            me%cont(j+nrd) = ydiff
                            bspl = h*k1(i) - ydiff
                            me%cont(j+nrd*2) = bspl
                            me%cont(j+nrd*3) = ydiff - h*k4(i) - bspl
                            me%cont(j+nrd*4) = d41*k1(i)+d46*k6(i)+d47*k7(i)+&
                                               d48*k8(i)+d49*k9(i)+d410*k10(i)+&
                                               d411*k2(i)+d412*k3(i)
                            me%cont(j+nrd*5) = d51*k1(i)+d56*k6(i)+d57*k7(i)+&
                                               d58*k8(i)+d59*k9(i)+d510*k10(i)+&
                                               d511*k2(i)+d512*k3(i)
                            me%cont(j+nrd*6) = d61*k1(i)+d66*k6(i)+d67*k7(i)+&
                                               d68*k8(i)+d69*k9(i)+d610*k10(i)+&
                                               d611*k2(i)+d612*k3(i)
                            me%cont(j+nrd*7) = d71*k1(i)+d76*k6(i)+d77*k7(i)+&
                                               d78*k8(i)+d79*k9(i)+d710*k10(i)+&
                                               d711*k2(i)+d712*k3(i)
                        end do
                        ! the next three function evaluations
                        y1 = y + h*(a141*k1+a147*k7+a148*k8+a149*k9+&
                                    a1410*k10+a1411*k2+a1412*k3+a1413*k4)
                        call me%fcn(x+c14*h,y1,k10)
                        y1 = y + h*(a151*k1+a156*k6+a157*k7+a158*k8+&
                                    a1511*k2+a1512*k3+a1513*k4+a1514*k10)
                        call me%fcn(x+c15*h,y1,k2)
                        y1 = y + h*(a161*k1+a166*k6+a167*k7+a168*k8+a169*k9+&
                                    a1613*k4+a1614*k10+a1615*k2)
                        call me%fcn(x+c16*h,y1,k3)
                        nfcn = nfcn + 3
                        ! final preparation
                        do j = 1 , nrd
                            i = me%icomp(j)
                            me%cont(j+nrd*4) = h*(me%cont(j+nrd*4)+d413*k4(i)+&
                                               d414*k10(i)+d415*k2(i)+d416*k3(i))
                            me%cont(j+nrd*5) = h*(me%cont(j+nrd*5)+d513*k4(i)+&
                                               d514*k10(i)+d515*k2(i)+d516*k3(i))
                            me%cont(j+nrd*6) = h*(me%cont(j+nrd*6)+d613*k4(i)+&
                                               d614*k10(i)+d615*k2(i)+d616*k3(i))
                            me%cont(j+nrd*7) = h*(me%cont(j+nrd*7)+d713*k4(i)+&
                                               d714*k10(i)+d715*k2(i)+d716*k3(i))
                        end do
                        me%hout = h
                    end if
                    k1 = k4
                    y = k5
                    me%xold = x
                    x = xph
                    if ( iout==1 .or. iout==2 .or. event ) then
                        call me%solout(naccpt+1,me%xold,x,y,irtrn,xout)
                        if ( irtrn<0 ) exit !abort
                    end if
                    ! normal exit
                    if ( last ) then
                        h = hnew
                        idid = 1
                        return
                    end if
                    if ( abs(hnew)>hmax ) hnew = posneg*hmax
                    if ( reject ) hnew = posneg*min(abs(hnew),abs(h))
                    reject = .false.
                else
                    ! step is rejected
                    hnew = h/min(facc1,fac11/safe)
                    reject = .true.
                    if ( naccpt>=1 ) nrejct = nrejct + 1
                    last = .false.
                end if
                h = hnew
            end if
        end do
    end if

    if ( iprint/=0 ) write (iprint,'(A,E18.4)') ' exit of dop853 at x=', x
    idid = 2

    end subroutine dp86co