SUBROUTINE rkfn78( N, FCN, X, Y, YP, XEND, EPS, HMAX, H )
C ----------------------------------------------------------------------
C
C     RKFN78 ( N, FCN, X, Y, YP, XEND, EPS, HMAX, H )
C
C     Numerical solution of a system of second order
C     ordinary differential equations  y"=f(x,y,dy)
C     This is an embedded nystroem method of order 7(8)  
C     due to Fehlberg with stepsize control
C
C     input parameters
C     ----------------
C     N           dimension of the system (N.LE.51)
C     FCN         name (external) of subroutine computing the
C                 second derivative F(X,Y,DY):
C                    SUBROUTINE FCN(X,Y,DY,F)
C                    REAL*8 X,Y(N),DY(N),F(N)
C                    F(1)=...   ETC.
C     X           initial x-value
C     XEND        final x-value (XEND.GT.X)
C     Y(N)        initial values for y
C     YP(N)       initial values for y'
C     EPS         local tolerance
C     HMAX        maximal stepsize
C     H           initial stepsize guess
C
C     output parameters
C     -----------------
C     Y(N)        solution at xend
C     YP(N)       derivative of solution at xend
C
C     COMMON STAT can be used for statistics
C        NFCN      number of function evaluations
C        NSTEP     number of computed steps
C        NACCPT    number of accepted steps
C        NREJCT    number of rejected steps
C
C  Ref.: Erwin Fehlberg, Computing 14, p.371, 1975
C ----------------------------------------------------------------------
C23456789012345678901234567890123456789012345678901234567890123456789012

      implicit none 

      ! I/O list

      INTEGER*4   N
      REAL*8      X, Y(N), YP(N), XEND, EPS, HMAX, H

      ! method related constants

      INTEGER*4   N_STAGE, P
      PARAMETER ( N_STAGE=13, P=7 )

      ! implementation dependant constants

      INTEGER*4   N_MAX, MAX_STEPS
      REAL*8      UROUND
      PARAMETER ( N_MAX=180, MAX_STEPS=2147483647, UROUND=1.1D-16 )
      ! auxiliary variables 

      REAL*8    K ( 1:N_MAX, 0:N_STAGE )
      REAL*8    Y1 ( 1:N_MAX ), YP1 ( 1:N_MAX ), TE ( 1:N_MAX )

      REAL*8    POSNEG, HNEW, DENOM, ERR, FAC, H2, SUM, SUMP
      INTEGER   I_STAGE, I, J, L 
      LOGICAL   REJECT

      ! coefficients

      REAL*8 ALPHA_(1:13), BETA_(1:13,0:12), GAMMA_(1:13,0:12)
      COMMON /CRKFN78/ ALPHA_, BETA_, GAMMA_ 


      ! statistics 

      INTEGER*4   NFCN,NSTEP,NACCPT,NREJCT
      COMMON    /STAT/ NFCN,NSTEP,NACCPT,NREJCT

Cf2py intent(in) N, FCN, X, XEND
Cf2py intent(inout) Y, YP, HMAX, H, EPS
Cf2py depend(in) Y, YP


      ! initial preparations

      IF (ALPHA_(13).NE.1.0D0) THEN
        CALL FN78INI
      END IF 

      POSNEG = SIGN ( 1.0D0, XEND-X)

      HMAX   = ABS(HMAX)
      H      = MIN ( MAX(1.0D-8,ABS(H)) , HMAX )
      H      = SIGN ( H, POSNEG )

      EPS    = MAX ( EPS, 9.0*UROUND )

      REJECT = .FALSE.

      NFCN   =  NFCN + 1

      CALL FCN ( X, Y, YP, K(1,0) )


      ! basic integration step 


      DO WHILE ( POSNEG*(X-XEND)+UROUND .LE. 0.0 )

        ! failure exit 

        IF ( NSTEP.GT.MAX_STEPS .OR. X+0.05*H.EQ.X ) THEN 
          WRITE (*,*) 'Exit of RKNF78 at x = ', X
          WRITE (*,*) 'NSTEP     = ', NSTEP
          WRITE (*,*) 'MAX_STEPS = ', MAX_STEPS
          WRITE (*,*) 'H = ', H
          WRITE (*,*) 'X+0.05*H = ', X+0.05*H
          STOP
        END IF

        ! limit step size

        IF ( (X+H-XEND)*POSNEG .GT. 0.0 ) THEN 
          H = XEND-X
        END IF
        H2    = H**2
        NSTEP = NSTEP+1

        ! calculate K(*,1) .. K(*,N_STAGE)

        DO I_STAGE = 1,N_STAGE 
          DO L=1,N
            SUM  = 0.0
            SUMP = 0.0
            DO J=0,(I_STAGE-1) 
              SUM  = SUM  + GAMMA_(I_STAGE,J)*K(L,J)
              SUMP = SUMP + BETA_ (I_STAGE,J)*K(L,J)
            END DO
            Y1(L)  = Y(L) + ALPHA_(I_STAGE)*H*YP(L) + H2*SUM
            YP1(L) = YP(L) + H*SUMP
          END DO
          CALL FCN( X + ALPHA_(I_STAGE)*H, Y1, YP1, K(1,I_STAGE) )
        END DO  ! end of loop over stages
        NFCN = NFCN+N_STAGE

        ! error term and relative error estimation

        DO L=1,N
          TE(L) = GAMMA_(N_STAGE,N_STAGE-1) * H2 
     .                  * ( K(L,N_STAGE-1) - K(L,N_STAGE) )
        END DO 
        ERR = 0.0
        DO L=1,N
          DENOM = MAX ( 1.0D-6, ABS(Y(L)), ABS(Y1(L)), 2.0*UROUND/EPS )
          ERR = ERR + ( TE(L) / DENOM )**2
        END DO
        ERR = SQRT ( ERR / N )

        ! new step size 0.1<HNEW/H<4.0 

        FAC = MAX( 0.25D0, MIN(10.0D0,(ERR/EPS)**(1.0/(P+2))/0.75D0) )
        HNEW = H / FAC

        ! acceptance test

        IF ( ERR .LT. EPS ) THEN
          ! step is accpeted;
          ! new solution 
          DO L=1,N
            SUM  = 0.0
            SUMP = 0.0
            DO I=0,N_STAGE-1   ! to be adapted 
              SUM  = SUM  + GAMMA_(N_STAGE,I)*K(L,I)
              SUMP = SUMP + BETA_(N_STAGE,I)*K(L,I)
            END DO 
            Y1(L)  = Y(L)  + H*YP(L) + H2*SUM    
            YP1(L) = YP(L) + H*SUMP
          END DO
          ! replace old values 
          X = X+H
          DO L=1,N
            Y(L)   = Y1(L)
            YP(L)  = YP1(L)
            K(L,0) = K(L,N_STAGE)
          END DO
          ! limit step size
          IF ( ABS(HNEW) .GT. HMAX ) THEN
            HNEW=POSNEG*HMAX
          END IF
          ! if last step was rejected then do not increase step size
          IF (REJECT) THEN
            HNEW=POSNEG*MIN(ABS(HNEW),ABS(H))
          END IF
          REJECT = .FALSE.
          NACCPT = NACCPT+1
        ELSE
          ! step is rejected
          REJECT = .TRUE.
          NREJCT = NREJCT+1                   
        ENDIF

        H = HNEW

      END DO  ! end of DO WHILE

      RETURN

      ! end of FILG11

      END


      !----------------------------------------------------------------
      !                                                                
      ! Subroutine FN78INI                                             
      !                                                                
      ! initialize constants for use in the runge-kuttta-nystrom       
      ! intgration routines RKNF78
      !                                                                
      !----------------------------------------------------------------

      SUBROUTINE FN78INI

      implicit none 

      REAL*8 ALPHA_(1:13), BETA_(1:13,0:12), GAMMA_(1:13,0:12)
      COMMON /CRKFN78/ ALPHA_, BETA_, GAMMA_ 

      ALPHA_(1) = 0.73470804841064383606103566874183D-1
      ALPHA_(2) = 0.11020620726159657540915525031127D0
      ALPHA_(3) = 0.16530931089239486311373302546691D0
      ALPHA_(4) = 0.5D0
      ALPHA_(5) = 0.26628826929126164263520369439706D0
      ALPHA_(6) = 0.63371173070873835736479630560294D0
      ALPHA_(7) = 0.75D0
      ALPHA_(8) = 0.5625D0
      ALPHA_(9) = 0.125D0
      ALPHA_(10) = 0.375D0
      ALPHA_(11) = 0.96652805160235700834451642119099D0
      ALPHA_(12) = 1.D0
      ALPHA_(13) = 1.D0

      BETA_(1,0) = 0.73470804841064383606103566874183D-1
      BETA_(2,0) = 0.27551551815399143852288837577819D-1
      BETA_(2,1) = 0.82654655446197431556866512733456D-1
      BETA_(3,0) = 0.41327327723098715778433256366728D-1
      BETA_(3,1) = 0.D0
      BETA_(3,2) = 0.12398198316929614733529976910018D0
      BETA_(4,0) = 0.89670558763795782658378325389875D0
      BETA_(4,1) = 0.D0
      BETA_(4,2) = -0.34585915268314336799411093327799D+1
      BETA_(4,3) = 0.30618859391934758533573260788811D+1
      BETA_(5,0) = 0.57053369423965328229363646644660D-1
      BETA_(5,1) = 0.D0
      BETA_(5,2) = 0.D0
      BETA_(5,3) = 0.20664670695498245793175075857561D0
      BETA_(5,4) = 0.25881929123138564740892891767848D-2
      BETA_(6,0) = 0.22130953402732738534298054285043D-1
      BETA_(6,1) = 0.D0
      BETA_(6,2) = 0.D0
      BETA_(6,3) = 0.36666901842159380713935315648106D0
      BETA_(6,4) = 0.32075560532767078702498038199180D0
      BETA_(6,5) = -0.75843846443258975333835287154967D-1
      BETA_(7,0) = 0.833333333333333333333333333333333D-1
      BETA_(7,1) = 0.D0
      BETA_(7,2) = 0.D0
      BETA_(7,3) = 0.D0
      BETA_(7,4) = 0.D0
      BETA_(7,5) = 0.38436436964131621037911008488971D0
      BETA_(7,6) = 0.28230229702535045628755658177969D0
      BETA_(8,0) = 0.8349609375D-1
      BETA_(8,1) = 0.D0
      BETA_(8,2) = 0.D0
      BETA_(8,3) = 0.D0
      BETA_(8,4) = 0.D0
      BETA_(8,5) = 0.38306747479397408845870063561136D0
      BETA_(8,6) = 0.13548721270602591154129936438864D0
      BETA_(8,7) = -0.3955078125D-1
      BETA_(9,0) = 0.73420353223593964334705075445816D-1
      BETA_(9,1) = 0.D0
      BETA_(9,2) = 0.D0
      BETA_(9,3) = 0.D0
      BETA_(9,4) = 0.D0
      BETA_(9,5) = 0.98808964916022916024205710336420D-1
      BETA_(9,6) = 0.24153311327327749549842803451955D0
      BETA_(9,7) = -0.48707561728395061728395061728395D-1
      BETA_(9,8) = -0.24005486968449931412894375857339D0
      BETA_(10,0) = 0.81378441127067064904041056404207D-2
      BETA_(10,1) = 0.D0
      BETA_(10,2) = 0.D0
      BETA_(10,3) = 0.D0
      BETA_(10,4) = 0.D0
      BETA_(10,5) = 0.D0
      BETA_(10,6) = -0.36266091174647134384031532058792D0
      BETA_(10,7) = 0.69726880597127928317272609847243D-1
      BETA_(10,8) = 0.37797780620763392161154341509711D0
      BETA_(10,9) = 0.28181838082900278742109519000315D0
      BETA_(11,0) = -0.14042538922482838913280031225476D+1
      BETA_(11,1) = 0.D0
      BETA_(11,2) = 0.D0
      BETA_(11,3) = 0.D0
      BETA_(11,4) = 0.D0
      BETA_(11,5) = -0.13555559029404957528304113342361D+2
      BETA_(11,6) = -0.15021472824848050961721330969968D+1
      BETA_(11,7) = 0.14767543284167949686233606841588D+1
      BETA_(11,8) = -0.21707681965133688432577373607995D+1
      BETA_(11,9) = 0.66149759502676558681039202833030D+1
      BETA_(11,10) = 0.11507526173569321530679222376434D+2
      BETA_(12,0) = -0.52708651815801315268176882187497D+1
      BETA_(12,1) = 0.D0
      BETA_(12,2) = 0.D0
      BETA_(12,3) = 0.D0
      BETA_(12,4) = 0.D0
      BETA_(12,5) = -0.49965599553656833001045921529105D+2
      BETA_(12,6) = -0.50302228928658231516135124812231D+1
      BETA_(12,7) = 0.44548269045298760506518238622704D+1
      BETA_(12,8) = -0.86071533124033841312406742989148D+1
      BETA_(12,9) = 0.23840410046372287590078676456468D+2
      BETA_(12,10) = 0.41711581466028388124069667164840D+2
      BETA_(12,11) = -0.13297747642437995408237095558512D0
      BETA_(13,0) = 0.35099303056581883152660173681744D-1
      BETA_(13,1) = 0.D0
      BETA_(13,2) = 0.D0
      BETA_(13,3) = 0.D0
      BETA_(13,4) = 0.D0
      BETA_(13,5) = 0.D0
      BETA_(13,6) = 0.D0
      BETA_(13,7) = 0.25223475276631606400638853417712D0
      BETA_(13,8) = 0.11840033306876549234162515364336D0
      BETA_(13,9) = 0.20258133611250929893187899871888D0
      BETA_(13,10) = 0.26757025259420140796393329272621D0
      BETA_(13,11) = 0.16586384510629873791268098150965D0
      BETA_(13,12) = -0.41749822704672884309167134456960D-1

      GAMMA_(1,0) = 0.26989795819968848329994970508715D-2
      GAMMA_(2,0) = 0.30363520297464954371244341822304D-2
      GAMMA_(2,1) = 0.30363520297464954371244341822304D-2
      GAMMA_(3,0) = 0.68317920669296147335299769100184D-2
      GAMMA_(3,1) = 0.D0
      GAMMA_(3,2) = 0.68317920669296147335299769100184D-2
      GAMMA_(4,0) = -0.10263757731977888994310872824217D-2
      GAMMA_(4,1) = 0.D0
      GAMMA_(4,2) = 0.D0
      GAMMA_(4,3) = 0.12602637577319778889943108728242D0
      GAMMA_(5,0) = 0.98909903843107417913313499064241D-2
      GAMMA_(5,1) = 0.D0
      GAMMA_(5,2) = 0.20401758759111349514170518498571D-1
      GAMMA_(5,3) = 0.50265147713328703261825104735338D-2
      GAMMA_(5,4) = 0.13545726631277755415728360014730D-3
      GAMMA_(6,0) = 0.36772464695317721429741572462201D-1
      GAMMA_(6,1) = 0.D0
      GAMMA_(6,2) = 0.D0
      GAMMA_(6,3) = 0.82132294778521785827721741407693D-1
      GAMMA_(6,4) = 0.30087165409098963036870918119641D-1
      GAMMA_(6,5) = 0.51803353935993790519824105531789D-1
      GAMMA_(7,0) = 0.41233049088272873123221004021091D-1
      GAMMA_(7,1) = 0.D0
      GAMMA_(7,2) = 0.D0
      GAMMA_(7,3) = 0.11335100293061819105328798078376D0
      GAMMA_(7,4) = 0.56722148592237668841301677436715D-1
      GAMMA_(7,5) = 0.57456202064954525469376924736474D-1
      GAMMA_(7,6) = 0.12487597323916741512812413021961D-1
      GAMMA_(8,0) = 0.4214630126953125D-1
      GAMMA_(8,1) = 0.D0
      GAMMA_(8,2) = 0.D0
      GAMMA_(8,3) = 0.D0
      GAMMA_(8,4) = -0.7808807373046875D-1
      GAMMA_(8,5) = 0.14104682102928772004397085536135D0
      GAMMA_(8,6) = 0.74603813736337279956029144638648D-1
      GAMMA_(8,7) = -0.215057373046875D-1
      GAMMA_(9,0) = 0.55243877171925011431184270690444D-2
      GAMMA_(9,1) = 0.D0
      GAMMA_(9,2) = 0.D0
      GAMMA_(9,3) = 0.D0
      GAMMA_(9,4) = 0.D0
      GAMMA_(9,5) = 0.45913375893505158838018111807029D-2
      GAMMA_(9,6) = 0.12009956992268139808927955623138D-1
      GAMMA_(9,7) = -0.24361818415637860082304526748971D-2
      GAMMA_(9,8) = -0.11877000457247370827617741197988D-1
      GAMMA_(10,0) = 0.12396099092300428073855581211091D-1
      GAMMA_(10,1) = 0.D0
      GAMMA_(10,2) = 0.D0
      GAMMA_(10,3) = 0.D0
      GAMMA_(10,4) = 0.D0
      GAMMA_(10,5) = 0.D0
      GAMMA_(10,6) = -0.23148568834881149606828637484335D-1
      GAMMA_(10,7) = 0.44057716338592294670599538200368D-2
      GAMMA_(10,8) = 0.24164236870396086181891828927171D-1
      GAMMA_(10,9) = 0.52494961238325405884021273526037D-1
      GAMMA_(11,0) = -0.12148292337172366838692371706654D0
      GAMMA_(11,1) = 0.D0
      GAMMA_(11,2) = 0.D0
      GAMMA_(11,3) = -0.15948786809469047245658868763595D+1
      GAMMA_(11,4) = 0.77089844409590354601143580630038D-1
      GAMMA_(11,5) = 0.D0
      GAMMA_(11,6) = 0.D0
      GAMMA_(11,7) = 0.98844932135442618048624328441900D-1
      GAMMA_(11,8) = -0.18517690177654009760124559303975D0
      GAMMA_(11,9) = 0.16665727117807342381867154630279D+1
      GAMMA_(11,10) = 0.52611925503652552568040599022802D0
      GAMMA_(12,0) = -0.49475846764102332689603709547782D0
      GAMMA_(12,1) = 0.D0
      GAMMA_(12,2) = 0.D0
      GAMMA_(12,3) = -0.56513209641364305307648232070852D+1
      GAMMA_(12,4) = 0.42750028729043677987389324310306D0
      GAMMA_(12,5) = 0.D0
      GAMMA_(12,6) = 0.D0
      GAMMA_(12,7) = 0.30293416726956828108567954376506D0
      GAMMA_(12,8) = -0.10280329379503342611614151091571D+1
      GAMMA_(12,9) = 0.54254171279669182157854764162220D+1
      GAMMA_(12,10) = 0.15340242607867031086671199895920D+1
      GAMMA_(12,11) = -0.15763473585838266589893780962028D-1
      GAMMA_(13,0) = 0.35173987586306713954725819037907D-1
      GAMMA_(13,1) = 0.D0
      GAMMA_(13,2) = 0.D0
      GAMMA_(13,3) = 0.D0
      GAMMA_(13,4) = 0.D0
      GAMMA_(13,5) = 0.D0
      GAMMA_(13,6) = 0.D0
      GAMMA_(13,7) = 0.63858784354258308506882892800822D-1
      GAMMA_(13,8) = 0.50866724905581448754291007148603D-1
      GAMMA_(13,9) = 0.17703179472766752427031494226269D0
      GAMMA_(13,10) = 0.16781715613041509463911067215393D0
      GAMMA_(13,11) = 0.45385629257942440722375392950401D-2
      GAMMA_(13,12) = 0.71298936997666580243712730101115D-3

      RETURN
      END