************************************************************************
************************************************************************
************************************************************************
      SUBROUTINE MMM(XIN,NVAR,YIN,DYDX,HTOT,NSTEP,DERIVS, YOUT)
C                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  ^^^^
C                                      in                   out


C     implements the modified midpoint method for use in Bulirsch-Stoer
C     integration of differential equations.

C     Input at XIN are NV dependent variables YIN and their derivatives
C     DYDX as well as the total step to be taken, HTOT, and the number
C     of substeps NSTEP.  DERIVS is an external routine for evaluating
C     the right hand sides.
C
C     The output is returned as YOUT.
C
C     Based on ideas from NUMERICAL RECIPES.

C     Code written by Phillip Helbig, Hamburg Observatory, Sep 1995.


************************************************************************
************************************************************************
C     declarations


************************************************************************
C     named constants for array bounds


      INTEGER NMAX

      PARAMETER (NMAX=50)



************************************************************************
C     input variables in the parameter list


C     dimension of vector and number of substeps to be used

      INTEGER NVAR, NSTEP

C     independent variable, input vector and deriviative vector, total
C     step to be taken

      DOUBLE PRECISION XIN, YIN(NVAR), DYDX(NVAR), HTOT



************************************************************************
C     output variables in the parameter list


C     output vector

      DOUBLE PRECISION YOUT(NVAR)



************************************************************************
C     output to called routines


C     all called routines

      EXTERNAL DERIVS

C     variables

C     x and y for external evaluation routine

      DOUBLE PRECISION X, Y2(NMAX)



************************************************************************
C     internally used quantities


C     variables

C     for DO loops

      INTEGER I, N

C     increment, 2 temporaries, second vector variable (use see below)

      DOUBLE PRECISION H, H2, SWAP, Y1(NMAX)



************************************************************************
************************************************************************
C     calculations


C     set the stepsize

      H = HTOT/NSTEP

C     initialise the work vectors for the first step

      DO 100 I = 1, NVAR

         Y1(I) = YIN(I)

         Y2(I) = YIN(I) + H*DYDX(I)

100   CONTINUE

C     increment X

      X = XIN + H

C     temporarily store derivatives in YOUT

      CALL DERIVS(X,Y2, YOUT)

C     just to save some work in the DO loop

      H2 = 2.0D0*H

C     basic step

      DO 300 N = 2, NSTEP

         DO 200 I = 1, NVAR

            SWAP = Y1(I) + H2*YOUT(I)

            Y1(I) = Y2(I)

            Y2(I) = SWAP

200      CONTINUE

C        increment X

         X = X + H

C        derivative evaluation

         CALL DERIVS(X,Y2, YOUT)

300   CONTINUE

C     this is the last step

      DO 400 I = 1, NVAR

         YOUT(I) = 0.5D0*(Y1(I) + Y2(I) + H*YOUT(I))

400   CONTINUE



************************************************************************
************************************************************************
      END