[Scipy-svn] r5589 - in trunk/scipy/integrate: . dop tests

scipy-svn@scip... scipy-svn@scip...
Mon Feb 23 04:26:46 CST 2009


Author: jtravs
Date: 2009-02-23 04:25:08 -0600 (Mon, 23 Feb 2009)
New Revision: 5589

Added:
   trunk/scipy/integrate/dop.pyf
   trunk/scipy/integrate/dop/
   trunk/scipy/integrate/dop/dop853.f
   trunk/scipy/integrate/dop/dopri5.f
Modified:
   trunk/scipy/integrate/ode.py
   trunk/scipy/integrate/setup.py
   trunk/scipy/integrate/tests/test_integrate.py
Log:
Add dopir5 and dop853 solvers to integrate.ode


Added: trunk/scipy/integrate/dop/dop853.f
===================================================================
--- trunk/scipy/integrate/dop/dop853.f	2009-02-23 05:02:34 UTC (rev 5588)
+++ trunk/scipy/integrate/dop/dop853.f	2009-02-23 10:25:08 UTC (rev 5589)
@@ -0,0 +1,879 @@
+      SUBROUTINE DOP853(N,FCN,X,Y,XEND,
+     &                  RTOL,ATOL,ITOL,
+     &                  SOLOUT,IOUT,
+     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
+C ----------------------------------------------------------
+C     NUMERICAL SOLUTION OF A SYSTEM OF FIRST 0RDER
+C     ORDINARY DIFFERENTIAL EQUATIONS  Y'=F(X,Y).
+C     THIS IS AN EXPLICIT RUNGE-KUTTA METHOD OF ORDER 8(5,3)  
+C     DUE TO DORMAND & PRINCE (WITH STEPSIZE CONTROL AND
+C     DENSE OUTPUT)
+C
+C     AUTHORS: E. HAIRER AND G. WANNER
+C              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
+C              CH-1211 GENEVE 24, SWITZERLAND 
+C              E-MAIL:  Ernst.Hairer@math.unige.ch
+C                       Gerhard.Wanner@math.unige.ch
+C     
+C     THIS CODE IS DESCRIBED IN:
+C         E. HAIRER, S.P. NORSETT AND G. WANNER, SOLVING ORDINARY
+C         DIFFERENTIAL EQUATIONS I. NONSTIFF PROBLEMS. 2ND EDITION. 
+C         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, 
+C         SPRINGER-VERLAG (1993)
+C      
+C     VERSION OF APRIL 25, 1996
+C     (latest correction of a small bug: August 8, 2005)
+C
+C     Edited (22 Feb 2009) by J.C. Travers:
+C       renamed HINIT->HINIT853 to avoid name collision with dopri5
+C
+C     INPUT PARAMETERS  
+C     ----------------  
+C     N           DIMENSION OF THE SYSTEM 
+C
+C     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
+C                 VALUE OF F(X,Y):
+C                    SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR)
+C                    DOUBLE PRECISION X,Y(N),F(N)
+C                    F(1)=...   ETC.
+C
+C     X           INITIAL X-VALUE
+C
+C     Y(N)        INITIAL VALUES FOR Y
+C
+C     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
+C
+C     RTOL,ATOL   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
+C                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
+C                 ATOL SHOULD BE STRICTLY POSITIVE (POSSIBLY VERY SMALL)
+C
+C     ITOL        SWITCH FOR RTOL AND ATOL:
+C                   ITOL=0: BOTH RTOL AND ATOL ARE SCALARS.
+C                     THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF
+C                     Y(I) BELOW RTOL*ABS(Y(I))+ATOL
+C                   ITOL=1: BOTH RTOL AND ATOL ARE VECTORS.
+C                     THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW
+C                     RTOL(I)*ABS(Y(I))+ATOL(I).
+C
+C     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
+C                 NUMERICAL SOLUTION DURING INTEGRATION. 
+C                 IF IOUT.GE.1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
+C                 SUPPLY A DUMMY SUBROUTINE IF IOUT=0. 
+C                 IT MUST HAVE THE FORM
+C                    SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,CON,ICOMP,ND,
+C                                       RPAR,IPAR,IRTRN)
+C                    DIMENSION Y(N),CON(8*ND),ICOMP(ND)
+C                    ....  
+C                 SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
+C                    GRID-POINT "X" (THEREBY THE INITIAL VALUE IS
+C                    THE FIRST GRID-POINT).
+C                 "XOLD" IS THE PRECEEDING GRID-POINT.
+C                 "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
+C                    IS SET <0, DOP853 WILL RETURN TO THE CALLING PROGRAM.
+C                    IF THE NUMERICAL SOLUTION IS ALTERED IN SOLOUT,
+C                    SET  IRTRN = 2
+C           
+C          -----  CONTINUOUS OUTPUT: -----
+C                 DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
+C                 FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH
+C                 THE FUNCTION
+C                        >>>   CONTD8(I,S,CON,ICOMP,ND)   <<<
+C                 WHICH PROVIDES AN APPROXIMATION TO THE I-TH
+C                 COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
+C                 S SHOULD LIE IN THE INTERVAL [XOLD,X].
+C           
+C     IOUT        SWITCH FOR CALLING THE SUBROUTINE SOLOUT:
+C                    IOUT=0: SUBROUTINE IS NEVER CALLED
+C                    IOUT=1: SUBROUTINE IS USED FOR OUTPUT
+C                    IOUT=2: DENSE OUTPUT IS PERFORMED IN SOLOUT
+C                            (IN THIS CASE WORK(5) MUST BE SPECIFIED)
+C
+C     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
+C                 WORK(1),...,WORK(20) SERVE AS PARAMETERS FOR THE CODE.
+C                 FOR STANDARD USE, SET THEM TO ZERO BEFORE CALLING.
+C                 "LWORK" MUST BE AT LEAST  11*N+8*NRDENS+21
+C                 WHERE  NRDENS = IWORK(5)
+C
+C     LWORK       DECLARED LENGHT OF ARRAY "WORK".
+C
+C     IWORK       INTEGER WORKING SPACE OF LENGHT "LIWORK".
+C                 IWORK(1),...,IWORK(20) SERVE AS PARAMETERS FOR THE CODE.
+C                 FOR STANDARD USE, SET THEM TO ZERO BEFORE CALLING.
+C                 "LIWORK" MUST BE AT LEAST NRDENS+21 .
+C
+C     LIWORK      DECLARED LENGHT OF ARRAY "IWORK".
+C
+C     RPAR, IPAR  REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH  
+C                 CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING
+C                 PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES. 
+C
+C-----------------------------------------------------------------------
+C 
+C     SOPHISTICATED SETTING OF PARAMETERS
+C     -----------------------------------
+C              SEVERAL PARAMETERS (WORK(1),...,IWORK(1),...) ALLOW
+C              TO ADAPT THE CODE TO THE PROBLEM AND TO THE NEEDS OF
+C              THE USER. FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES.
+C
+C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 2.3D-16.
+C
+C    WORK(2)   THE SAFETY FACTOR IN STEP SIZE PREDICTION,
+C              DEFAULT 0.9D0.
+C
+C    WORK(3), WORK(4)   PARAMETERS FOR STEP SIZE SELECTION
+C              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION
+C                 WORK(3) <= HNEW/HOLD <= WORK(4)
+C              DEFAULT VALUES: WORK(3)=0.333D0, WORK(4)=6.D0
+C
+C    WORK(5)   IS THE "BETA" FOR STABILIZED STEP SIZE CONTROL
+C              (SEE SECTION IV.2). POSITIVE VALUES OF BETA ( <= 0.04 )
+C              MAKE THE STEP SIZE CONTROL MORE STABLE.
+C              NEGATIVE WORK(5) PROVOKE BETA=0.
+C              DEFAULT 0.0D0.
+C
+C    WORK(6)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
+C
+C    WORK(7)   INITIAL STEP SIZE, FOR WORK(7)=0.D0 AN INITIAL GUESS
+C              IS COMPUTED WITH HELP OF THE FUNCTION HINIT
+C
+C    IWORK(1)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
+C              THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000.
+C
+C    IWORK(2)  SWITCH FOR THE CHOICE OF THE COEFFICIENTS
+C              IF IWORK(2).EQ.1  METHOD DOP853 OF DORMAND AND PRINCE
+C              (SECTION II.6).
+C              THE DEFAULT VALUE (FOR IWORK(2)=0) IS IWORK(2)=1.
+C
+C    IWORK(3)  SWITCH FOR PRINTING ERROR MESSAGES
+C              IF IWORK(3).LT.0 NO MESSAGES ARE BEING PRINTED
+C              IF IWORK(3).GT.0 MESSAGES ARE PRINTED WITH
+C              WRITE (IWORK(3),*) ...  
+C              DEFAULT VALUE (FOR IWORK(3)=0) IS IWORK(3)=6
+C
+C    IWORK(4)  TEST FOR STIFFNESS IS ACTIVATED AFTER STEP NUMBER
+C              J*IWORK(4) (J INTEGER), PROVIDED IWORK(4).GT.0.
+C              FOR NEGATIVE IWORK(4) THE STIFFNESS TEST IS
+C              NEVER ACTIVATED; DEFAULT VALUE IS IWORK(4)=1000
+C
+C    IWORK(5)  = NRDENS = NUMBER OF COMPONENTS, FOR WHICH DENSE OUTPUT
+C              IS REQUIRED; DEFAULT VALUE IS IWORK(5)=0;
+C              FOR   0 < NRDENS < N   THE COMPONENTS (FOR WHICH DENSE
+C              OUTPUT IS REQUIRED) HAVE TO BE SPECIFIED IN
+C              IWORK(21),...,IWORK(NRDENS+20);
+C              FOR  NRDENS=N  THIS IS DONE BY THE CODE.
+C
+C----------------------------------------------------------------------
+C
+C     OUTPUT PARAMETERS 
+C     ----------------- 
+C     X           X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED
+C                 (AFTER SUCCESSFUL RETURN X=XEND).
+C
+C     Y(N)        NUMERICAL SOLUTION AT X
+C 
+C     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
+C
+C     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
+C                   IDID= 1  COMPUTATION SUCCESSFUL,
+C                   IDID= 2  COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT)
+C                   IDID=-1  INPUT IS NOT CONSISTENT,
+C                   IDID=-2  LARGER NMAX IS NEEDED,
+C                   IDID=-3  STEP SIZE BECOMES TOO SMALL.
+C                   IDID=-4  PROBLEM IS PROBABLY STIFF (INTERRUPTED).
+C
+C   IWORK(17)  NFCN    NUMBER OF FUNCTION EVALUATIONS
+C   IWORK(18)  NSTEP   NUMBER OF COMPUTED STEPS
+C   IWORK(19)  NACCPT  NUMBER OF ACCEPTED STEPS
+C   IWORK(20)  NREJCT  NUMBER OF REJECTED STEPS (DUE TO ERROR TEST),
+C                      (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED)
+C-----------------------------------------------------------------------
+C *** *** *** *** *** *** *** *** *** *** *** *** ***
+C          DECLARATIONS 
+C *** *** *** *** *** *** *** *** *** *** *** *** ***
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      DIMENSION Y(N),ATOL(*),RTOL(*),WORK(LWORK),IWORK(LIWORK)
+      DIMENSION RPAR(*),IPAR(*)
+      LOGICAL ARRET
+      EXTERNAL FCN,SOLOUT
+C *** *** *** *** *** *** ***
+C        SETTING THE PARAMETERS 
+C *** *** *** *** *** *** ***
+      NFCN=0
+      NSTEP=0
+      NACCPT=0
+      NREJCT=0
+      ARRET=.FALSE.
+C -------- IPRINT FOR MONITORING THE PRINTING
+      IF(IWORK(3).EQ.0)THEN
+         IPRINT=6
+      ELSE
+         IPRINT=IWORK(3)
+      END IF
+C -------- NMAX , THE MAXIMAL NUMBER OF STEPS ----- 
+      IF(IWORK(1).EQ.0)THEN
+         NMAX=100000
+      ELSE
+         NMAX=IWORK(1)
+         IF(NMAX.LE.0)THEN
+            IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &          ' WRONG INPUT IWORK(1)=',IWORK(1)
+            ARRET=.TRUE.
+         END IF
+      END IF
+C -------- METH   COEFFICIENTS OF THE METHOD
+      IF(IWORK(2).EQ.0)THEN
+         METH=1
+      ELSE
+         METH=IWORK(2)
+         IF(METH.LE.0.OR.METH.GE.4)THEN
+            IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &          ' CURIOUS INPUT IWORK(2)=',IWORK(2)
+            ARRET=.TRUE.
+         END IF
+      END IF  
+C -------- NSTIFF   PARAMETER FOR STIFFNESS DETECTION  
+      NSTIFF=IWORK(4) 
+      IF (NSTIFF.EQ.0) NSTIFF=1000
+      IF (NSTIFF.LT.0) NSTIFF=NMAX+10
+C -------- NRDENS   NUMBER OF DENSE OUTPUT COMPONENTS
+      NRDENS=IWORK(5)
+      IF(NRDENS.LT.0.OR.NRDENS.GT.N)THEN
+         IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &           ' CURIOUS INPUT IWORK(5)=',IWORK(5)
+         ARRET=.TRUE.
+      ELSE
+         IF(NRDENS.GT.0.AND.IOUT.LT.2)THEN
+            IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &       ' WARNING: PUT IOUT=2 FOR DENSE OUTPUT '
+         END IF 
+         IF (NRDENS.EQ.N) THEN
+            DO I=1,NRDENS
+               IWORK(I+20)=I
+            END DO
+         END IF
+      END IF       
+C -------- UROUND   SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0  
+      IF(WORK(1).EQ.0.D0)THEN
+         UROUND=2.3D-16
+      ELSE
+         UROUND=WORK(1)
+         IF(UROUND.LE.1.D-35.OR.UROUND.GE.1.D0)THEN
+            IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &        ' WHICH MACHINE DO YOU HAVE? YOUR UROUND WAS:',WORK(1)
+            ARRET=.TRUE.
+         END IF
+      END IF
+C -------  SAFETY FACTOR -------------
+      IF(WORK(2).EQ.0.D0)THEN
+         SAFE=0.9D0
+      ELSE
+         SAFE=WORK(2)
+         IF(SAFE.GE.1.D0.OR.SAFE.LE.1.D-4)THEN
+            IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &          ' CURIOUS INPUT FOR SAFETY FACTOR WORK(2)=',WORK(2)
+            ARRET=.TRUE.
+         END IF
+      END IF
+C -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
+      IF(WORK(3).EQ.0.D0)THEN
+         FAC1=0.333D0
+      ELSE
+         FAC1=WORK(3)
+      END IF
+      IF(WORK(4).EQ.0.D0)THEN
+         FAC2=6.D0
+      ELSE
+         FAC2=WORK(4)
+      END IF
+C --------- BETA FOR STEP CONTROL STABILIZATION -----------
+      IF(WORK(5).EQ.0.D0)THEN
+         BETA=0.0D0
+      ELSE
+         IF(WORK(5).LT.0.D0)THEN
+            BETA=0.D0
+         ELSE
+            BETA=WORK(5)
+            IF(BETA.GT.0.2D0)THEN
+               IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &          ' CURIOUS INPUT FOR BETA: WORK(5)=',WORK(5)
+            ARRET=.TRUE.
+         END IF
+         END IF
+      END IF
+C -------- MAXIMAL STEP SIZE
+      IF(WORK(6).EQ.0.D0)THEN
+         HMAX=XEND-X
+      ELSE
+         HMAX=WORK(6)
+      END IF
+C -------- INITIAL STEP SIZE
+      H=WORK(7)
+C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
+      IEK1=21
+      IEK2=IEK1+N
+      IEK3=IEK2+N
+      IEK4=IEK3+N
+      IEK5=IEK4+N
+      IEK6=IEK5+N
+      IEK7=IEK6+N
+      IEK8=IEK7+N
+      IEK9=IEK8+N
+      IEK10=IEK9+N
+      IEY1=IEK10+N
+      IECO=IEY1+N
+C ------ TOTAL STORAGE REQUIREMENT -----------
+      ISTORE=IECO+8*NRDENS-1
+      IF(ISTORE.GT.LWORK)THEN
+        IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &   ' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
+        ARRET=.TRUE.
+      END IF
+      ICOMP=21
+      ISTORE=ICOMP+NRDENS-1
+      IF(ISTORE.GT.LIWORK)THEN
+        IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &   ' INSUFFICIENT STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
+        ARRET=.TRUE.
+      END IF
+C -------- WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
+      IF (ARRET) THEN
+         IDID=-1
+         RETURN
+      END IF
+C -------- CALL TO CORE INTEGRATOR ------------
+      CALL DP86CO(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,IPRINT,
+     &   SOLOUT,IOUT,IDID,NMAX,UROUND,METH,NSTIFF,SAFE,BETA,FAC1,FAC2,
+     &   WORK(IEK1),WORK(IEK2),WORK(IEK3),WORK(IEK4),WORK(IEK5),
+     &   WORK(IEK6),WORK(IEK7),WORK(IEK8),WORK(IEK9),WORK(IEK10),
+     &   WORK(IEY1),WORK(IECO),IWORK(ICOMP),NRDENS,RPAR,IPAR,
+     &   NFCN,NSTEP,NACCPT,NREJCT)
+      WORK(7)=H
+      IWORK(17)=NFCN
+      IWORK(18)=NSTEP
+      IWORK(19)=NACCPT
+      IWORK(20)=NREJCT
+C ----------- RETURN -----------
+      RETURN
+      END
+C
+C
+C
+C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
+C
+      SUBROUTINE DP86CO(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,IPRINT,
+     &   SOLOUT,IOUT,IDID,NMAX,UROUND,METH,NSTIFF,SAFE,BETA,FAC1,FAC2,
+     &   K1,K2,K3,K4,K5,K6,K7,K8,K9,K10,Y1,CONT,ICOMP,NRD,RPAR,IPAR,
+     &   NFCN,NSTEP,NACCPT,NREJCT)
+C ----------------------------------------------------------
+C     CORE INTEGRATOR FOR DOP853
+C     PARAMETERS SAME AS IN DOP853 WITH WORKSPACE ADDED 
+C ---------------------------------------------------------- 
+C         DECLARATIONS 
+C ---------------------------------------------------------- 
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      parameter (
+     &  c2  = 0.526001519587677318785587544488D-01,
+     &  c3  = 0.789002279381515978178381316732D-01,
+     &  c4  = 0.118350341907227396726757197510D+00,
+     &  c5  = 0.281649658092772603273242802490D+00,
+     &  c6  = 0.333333333333333333333333333333D+00,
+     &  c7  = 0.25D+00,
+     &  c8  = 0.307692307692307692307692307692D+00,
+     &  c9  = 0.651282051282051282051282051282D+00,
+     &  c10 = 0.6D+00,
+     &  c11 = 0.857142857142857142857142857142D+00,
+     &  c14 = 0.1D+00,
+     &  c15 = 0.2D+00,
+     &  c16 = 0.777777777777777777777777777778D+00)
+      parameter (
+     &  b1 =   5.42937341165687622380535766363D-2,
+     &  b6 =   4.45031289275240888144113950566D0,
+     &  b7 =   1.89151789931450038304281599044D0,
+     &  b8 =  -5.8012039600105847814672114227D0,
+     &  b9 =   3.1116436695781989440891606237D-1,
+     &  b10 = -1.52160949662516078556178806805D-1,
+     &  b11 =  2.01365400804030348374776537501D-1,
+     &  b12 =  4.47106157277725905176885569043D-2)
+      parameter (
+     &  bhh1 = 0.244094488188976377952755905512D+00,
+     &  bhh2 = 0.733846688281611857341361741547D+00,
+     &  bhh3 = 0.220588235294117647058823529412D-01)
+      parameter (
+     &  er 1 =  0.1312004499419488073250102996D-01,
+     &  er 6 = -0.1225156446376204440720569753D+01,
+     &  er 7 = -0.4957589496572501915214079952D+00,
+     &  er 8 =  0.1664377182454986536961530415D+01,
+     &  er 9 = -0.3503288487499736816886487290D+00,
+     &  er10 =  0.3341791187130174790297318841D+00,
+     &  er11 =  0.8192320648511571246570742613D-01,
+     &  er12 = -0.2235530786388629525884427845D-01)
+      parameter (
+     &  a21 =    5.26001519587677318785587544488D-2,
+     &  a31 =    1.97250569845378994544595329183D-2,
+     &  a32 =    5.91751709536136983633785987549D-2,
+     &  a41 =    2.95875854768068491816892993775D-2,
+     &  a43 =    8.87627564304205475450678981324D-2,
+     &  a51 =    2.41365134159266685502369798665D-1,
+     &  a53 =   -8.84549479328286085344864962717D-1,
+     &  a54 =    9.24834003261792003115737966543D-1,
+     &  a61 =    3.7037037037037037037037037037D-2,
+     &  a64 =    1.70828608729473871279604482173D-1,
+     &  a65 =    1.25467687566822425016691814123D-1,
+     &  a71 =    3.7109375D-2,
+     &  a74 =    1.70252211019544039314978060272D-1,
+     &  a75 =    6.02165389804559606850219397283D-2,
+     &  a76 =   -1.7578125D-2)
+      parameter (
+     &  a81 =    3.70920001185047927108779319836D-2,
+     &  a84 =    1.70383925712239993810214054705D-1,
+     &  a85 =    1.07262030446373284651809199168D-1,
+     &  a86 =   -1.53194377486244017527936158236D-2,
+     &  a87 =    8.27378916381402288758473766002D-3,
+     &  a91 =    6.24110958716075717114429577812D-1,
+     &  a94 =   -3.36089262944694129406857109825D0,
+     &  a95 =   -8.68219346841726006818189891453D-1,
+     &  a96 =    2.75920996994467083049415600797D1,
+     &  a97 =    2.01540675504778934086186788979D1,
+     &  a98 =   -4.34898841810699588477366255144D1,
+     &  a101 =   4.77662536438264365890433908527D-1,
+     &  a104 =  -2.48811461997166764192642586468D0,
+     &  a105 =  -5.90290826836842996371446475743D-1,
+     &  a106 =   2.12300514481811942347288949897D1,
+     &  a107 =   1.52792336328824235832596922938D1,
+     &  a108 =  -3.32882109689848629194453265587D1,
+     &  a109 =  -2.03312017085086261358222928593D-2)
+      parameter (
+     &  a111 =  -9.3714243008598732571704021658D-1,
+     &  a114 =   5.18637242884406370830023853209D0,
+     &  a115 =   1.09143734899672957818500254654D0,
+     &  a116 =  -8.14978701074692612513997267357D0,
+     &  a117 =  -1.85200656599969598641566180701D1,
+     &  a118 =   2.27394870993505042818970056734D1,
+     &  a119 =   2.49360555267965238987089396762D0,
+     &  a1110 = -3.0467644718982195003823669022D0,
+     &  a121 =   2.27331014751653820792359768449D0,
+     &  a124 =  -1.05344954667372501984066689879D1,
+     &  a125 =  -2.00087205822486249909675718444D0,
+     &  a126 =  -1.79589318631187989172765950534D1,
+     &  a127 =   2.79488845294199600508499808837D1,
+     &  a128 =  -2.85899827713502369474065508674D0,
+     &  a129 =  -8.87285693353062954433549289258D0,
+     &  a1210 =  1.23605671757943030647266201528D1,
+     &  a1211 =  6.43392746015763530355970484046D-1)
+      parameter (
+     &  a141 =  5.61675022830479523392909219681D-2,
+     &  a147 =  2.53500210216624811088794765333D-1,
+     &  a148 = -2.46239037470802489917441475441D-1,
+     &  a149 = -1.24191423263816360469010140626D-1,
+     &  a1410 =  1.5329179827876569731206322685D-1,
+     &  a1411 =  8.20105229563468988491666602057D-3,
+     &  a1412 =  7.56789766054569976138603589584D-3,
+     &  a1413 = -8.298D-3)
+      parameter (
+     &  a151 =  3.18346481635021405060768473261D-2,
+     &  a156 =  2.83009096723667755288322961402D-2,
+     &  a157 =  5.35419883074385676223797384372D-2,
+     &  a158 = -5.49237485713909884646569340306D-2,
+     &  a1511 = -1.08347328697249322858509316994D-4,
+     &  a1512 =  3.82571090835658412954920192323D-4,
+     &  a1513 = -3.40465008687404560802977114492D-4,
+     &  a1514 =  1.41312443674632500278074618366D-1,
+     &  a161 = -4.28896301583791923408573538692D-1,
+     &  a166 = -4.69762141536116384314449447206D0,
+     &  a167 =  7.68342119606259904184240953878D0,
+     &  a168 =  4.06898981839711007970213554331D0,
+     &  a169 =  3.56727187455281109270669543021D-1,
+     &  a1613 = -1.39902416515901462129418009734D-3,
+     &  a1614 =  2.9475147891527723389556272149D0,
+     &  a1615 = -9.15095847217987001081870187138D0)
+      parameter (
+     &  d41  = -0.84289382761090128651353491142D+01,
+     &  d46  =  0.56671495351937776962531783590D+00,
+     &  d47  = -0.30689499459498916912797304727D+01,
+     &  d48  =  0.23846676565120698287728149680D+01,
+     &  d49  =  0.21170345824450282767155149946D+01,
+     &  d410 = -0.87139158377797299206789907490D+00,
+     &  d411 =  0.22404374302607882758541771650D+01,
+     &  d412 =  0.63157877876946881815570249290D+00,
+     &  d413 = -0.88990336451333310820698117400D-01,
+     &  d414 =  0.18148505520854727256656404962D+02,
+     &  d415 = -0.91946323924783554000451984436D+01,
+     &  d416 = -0.44360363875948939664310572000D+01)
+      parameter (
+     &  d51  =  0.10427508642579134603413151009D+02,
+     &  d56  =  0.24228349177525818288430175319D+03,
+     &  d57  =  0.16520045171727028198505394887D+03,
+     &  d58  = -0.37454675472269020279518312152D+03,
+     &  d59  = -0.22113666853125306036270938578D+02,
+     &  d510 =  0.77334326684722638389603898808D+01,
+     &  d511 = -0.30674084731089398182061213626D+02,
+     &  d512 = -0.93321305264302278729567221706D+01,
+     &  d513 =  0.15697238121770843886131091075D+02,
+     &  d514 = -0.31139403219565177677282850411D+02,
+     &  d515 = -0.93529243588444783865713862664D+01,
+     &  d516 =  0.35816841486394083752465898540D+02)
+      parameter (
+     &  d61 =  0.19985053242002433820987653617D+02,
+     &  d66 = -0.38703730874935176555105901742D+03,
+     &  d67 = -0.18917813819516756882830838328D+03,
+     &  d68 =  0.52780815920542364900561016686D+03,
+     &  d69 = -0.11573902539959630126141871134D+02,
+     &  d610 =  0.68812326946963000169666922661D+01,
+     &  d611 = -0.10006050966910838403183860980D+01,
+     &  d612 =  0.77771377980534432092869265740D+00,
+     &  d613 = -0.27782057523535084065932004339D+01,
+     &  d614 = -0.60196695231264120758267380846D+02,
+     &  d615 =  0.84320405506677161018159903784D+02,
+     &  d616 =  0.11992291136182789328035130030D+02)
+      parameter (
+     &  d71  = -0.25693933462703749003312586129D+02,
+     &  d76  = -0.15418974869023643374053993627D+03,
+     &  d77  = -0.23152937917604549567536039109D+03,
+     &  d78  =  0.35763911791061412378285349910D+03,
+     &  d79  =  0.93405324183624310003907691704D+02,
+     &  d710 = -0.37458323136451633156875139351D+02,
+     &  d711 =  0.10409964950896230045147246184D+03,
+     &  d712 =  0.29840293426660503123344363579D+02,
+     &  d713 = -0.43533456590011143754432175058D+02,
+     &  d714 =  0.96324553959188282948394950600D+02,
+     &  d715 = -0.39177261675615439165231486172D+02,
+     &  d716 = -0.14972683625798562581422125276D+03)
+      DOUBLE PRECISION Y(N),Y1(N),K1(N),K2(N),K3(N),K4(N),K5(N),K6(N)
+      DOUBLE PRECISION K7(N),K8(N),K9(N),K10(N),ATOL(*),RTOL(*)     
+      DIMENSION CONT(8*NRD),ICOMP(NRD),RPAR(*),IPAR(*)
+      LOGICAL REJECT,LAST 
+      EXTERNAL FCN
+      COMMON /CONDO8/XOLD,HOUT
+C *** *** *** *** *** *** ***
+C  INITIALISATIONS
+C *** *** *** *** *** *** *** 
+      FACOLD=1.D-4  
+      EXPO1=1.d0/8.d0-BETA*0.2D0
+      FACC1=1.D0/FAC1
+      FACC2=1.D0/FAC2
+      POSNEG=SIGN(1.D0,XEND-X) 
+C --- INITIAL PREPARATIONS   
+      ATOLI=ATOL(1)
+      RTOLI=RTOL(1)    
+      LAST=.FALSE. 
+      HLAMB=0.D0
+      IASTI=0
+      CALL FCN(N,X,Y,K1,RPAR,IPAR)
+      HMAX=ABS(HMAX)     
+      IORD=8  
+      IF (H.EQ.0.D0) H=HINIT853(N,FCN,X,Y,XEND,POSNEG,K1,K2,K3,IORD,
+     &                       HMAX,ATOL,RTOL,ITOL,RPAR,IPAR)
+      NFCN=NFCN+2
+      REJECT=.FALSE.
+      XOLD=X
+      IF (IOUT.GE.1) THEN 
+          IRTRN=1 
+          HOUT=1.D0
+          CALL SOLOUT(NACCPT+1,XOLD,X,Y,N,CONT,ICOMP,NRD,
+     &                RPAR,IPAR,IRTRN)
+          IF (IRTRN.LT.0) GOTO 79
+      END IF
+C --- BASIC INTEGRATION STEP  
+   1  CONTINUE
+      IF (NSTEP.GT.NMAX) GOTO 78
+      IF (0.1D0*ABS(H).LE.ABS(X)*UROUND)GOTO 77
+      IF ((X+1.01D0*H-XEND)*POSNEG.GT.0.D0) THEN
+         H=XEND-X 
+         LAST=.TRUE.
+      END IF
+      NSTEP=NSTEP+1
+C --- THE TWELVE STAGES
+      IF (IRTRN.GE.2) THEN
+         CALL FCN(N,X,Y,K1,RPAR,IPAR)
+      END IF
+      DO 22 I=1,N 
+  22  Y1(I)=Y(I)+H*A21*K1(I)  
+      CALL FCN(N,X+C2*H,Y1,K2,RPAR,IPAR)
+      DO 23 I=1,N 
+  23  Y1(I)=Y(I)+H*(A31*K1(I)+A32*K2(I))  
+      CALL FCN(N,X+C3*H,Y1,K3,RPAR,IPAR)
+      DO 24 I=1,N 
+  24  Y1(I)=Y(I)+H*(A41*K1(I)+A43*K3(I))  
+      CALL FCN(N,X+C4*H,Y1,K4,RPAR,IPAR)
+      DO 25 I=1,N 
+  25  Y1(I)=Y(I)+H*(A51*K1(I)+A53*K3(I)+A54*K4(I))
+      CALL FCN(N,X+C5*H,Y1,K5,RPAR,IPAR)
+      DO 26 I=1,N 
+  26  Y1(I)=Y(I)+H*(A61*K1(I)+A64*K4(I)+A65*K5(I))
+      CALL FCN(N,X+C6*H,Y1,K6,RPAR,IPAR)
+      DO 27 I=1,N 
+  27  Y1(I)=Y(I)+H*(A71*K1(I)+A74*K4(I)+A75*K5(I)+A76*K6(I))
+      CALL FCN(N,X+C7*H,Y1,K7,RPAR,IPAR)
+      DO 28 I=1,N 
+  28  Y1(I)=Y(I)+H*(A81*K1(I)+A84*K4(I)+A85*K5(I)+A86*K6(I)+A87*K7(I))  
+      CALL FCN(N,X+C8*H,Y1,K8,RPAR,IPAR)
+      DO 29 I=1,N 
+  29  Y1(I)=Y(I)+H*(A91*K1(I)+A94*K4(I)+A95*K5(I)+A96*K6(I)+A97*K7(I)
+     &   +A98*K8(I))
+      CALL FCN(N,X+C9*H,Y1,K9,RPAR,IPAR)
+      DO 30 I=1,N 
+  30  Y1(I)=Y(I)+H*(A101*K1(I)+A104*K4(I)+A105*K5(I)+A106*K6(I)
+     &   +A107*K7(I)+A108*K8(I)+A109*K9(I))
+      CALL FCN(N,X+C10*H,Y1,K10,RPAR,IPAR)
+      DO 31 I=1,N 
+  31  Y1(I)=Y(I)+H*(A111*K1(I)+A114*K4(I)+A115*K5(I)+A116*K6(I)
+     &   +A117*K7(I)+A118*K8(I)+A119*K9(I)+A1110*K10(I))
+      CALL FCN(N,X+C11*H,Y1,K2,RPAR,IPAR)
+      XPH=X+H
+      DO 32 I=1,N 
+  32  Y1(I)=Y(I)+H*(A121*K1(I)+A124*K4(I)+A125*K5(I)+A126*K6(I)
+     &   +A127*K7(I)+A128*K8(I)+A129*K9(I)+A1210*K10(I)+A1211*K2(I))
+      CALL FCN(N,XPH,Y1,K3,RPAR,IPAR)
+      NFCN=NFCN+11
+      DO 35 I=1,N 
+      K4(I)=B1*K1(I)+B6*K6(I)+B7*K7(I)+B8*K8(I)+B9*K9(I)
+     &   +B10*K10(I)+B11*K2(I)+B12*K3(I)
+  35  K5(I)=Y(I)+H*K4(I)
+C --- ERROR ESTIMATION  
+      ERR=0.D0
+      ERR2=0.D0
+      IF (ITOL.EQ.0) THEN   
+        DO 41 I=1,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)
+  41    ERR=ERR+(ERRI/SK)**2
+      ELSE
+        DO 42 I=1,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)
+  42    ERR=ERR+(ERRI/SK)**2
+      END IF  
+      DENO=ERR+0.01D0*ERR2
+      IF (DENO.LE.0.D0) DENO=1.D0
+      ERR=ABS(H)*ERR*SQRT(1.D0/(N*DENO))
+C --- COMPUTATION OF HNEW
+      FAC11=ERR**EXPO1
+C --- LUND-STABILIZATION
+      FAC=FAC11/FACOLD**BETA
+C --- WE REQUIRE  FAC1 <= HNEW/H <= FAC2
+      FAC=MAX(FACC2,MIN(FACC1,FAC/SAFE))
+      HNEW=H/FAC  
+      IF(ERR.LE.1.D0)THEN
+C --- STEP IS ACCEPTED  
+         FACOLD=MAX(ERR,1.0D-4)
+         NACCPT=NACCPT+1
+         CALL FCN(N,XPH,K5,K4,RPAR,IPAR)
+         NFCN=NFCN+1
+C ------- STIFFNESS DETECTION
+         IF (MOD(NACCPT,NSTIFF).EQ.0.OR.IASTI.GT.0) THEN
+            STNUM=0.D0
+            STDEN=0.D0
+            DO 64 I=1,N 
+               STNUM=STNUM+(K4(I)-K3(I))**2
+               STDEN=STDEN+(K5(I)-Y1(I))**2
+ 64         CONTINUE  
+            IF (STDEN.GT.0.D0) HLAMB=ABS(H)*SQRT(STNUM/STDEN) 
+            IF (HLAMB.GT.6.1D0) THEN
+               NONSTI=0
+               IASTI=IASTI+1  
+               IF (IASTI.EQ.15) THEN
+                  IF (IPRINT.GT.0) WRITE (IPRINT,*) 
+     &               ' THE PROBLEM SEEMS TO BECOME STIFF AT X = ',X   
+                  IF (IPRINT.LE.0) GOTO 76
+               END IF
+            ELSE
+               NONSTI=NONSTI+1  
+               IF (NONSTI.EQ.6) IASTI=0
+            END IF
+         END IF 
+C ------- FINAL PREPARATION FOR DENSE OUTPUT
+         IF (IOUT.GE.2) THEN
+C ----    SAVE THE FIRST FUNCTION EVALUATIONS   
+            DO 62 J=1,NRD
+               I=ICOMP(J)
+               CONT(J)=Y(I)
+               YDIFF=K5(I)-Y(I)
+               CONT(J+NRD)=YDIFF
+               BSPL=H*K1(I)-YDIFF
+               CONT(J+NRD*2)=BSPL
+               CONT(J+NRD*3)=YDIFF-H*K4(I)-BSPL
+               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)
+               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)
+               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)
+               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)
+   62       CONTINUE 
+C ---     THE NEXT THREE FUNCTION EVALUATIONS
+            DO 51 I=1,N 
+  51           Y1(I)=Y(I)+H*(A141*K1(I)+A147*K7(I)+A148*K8(I)
+     &            +A149*K9(I)+A1410*K10(I)+A1411*K2(I)+A1412*K3(I)
+     &            +A1413*K4(I))
+            CALL FCN(N,X+C14*H,Y1,K10,RPAR,IPAR)
+            DO 52 I=1,N 
+  52           Y1(I)=Y(I)+H*(A151*K1(I)+A156*K6(I)+A157*K7(I)
+     &            +A158*K8(I)+A1511*K2(I)+A1512*K3(I)+A1513*K4(I)
+     &            +A1514*K10(I))
+            CALL FCN(N,X+C15*H,Y1,K2,RPAR,IPAR)
+            DO 53 I=1,N 
+  53           Y1(I)=Y(I)+H*(A161*K1(I)+A166*K6(I)+A167*K7(I)
+     &            +A168*K8(I)+A169*K9(I)+A1613*K4(I)+A1614*K10(I)
+     &            +A1615*K2(I))
+            CALL FCN(N,X+C16*H,Y1,K3,RPAR,IPAR)
+            NFCN=NFCN+3 
+C ---     FINAL PREPARATION
+            DO 63 J=1,NRD
+               I=ICOMP(J)
+               CONT(J+NRD*4)=H*(CONT(J+NRD*4)+D413*K4(I)+D414*K10(I)
+     &            +D415*K2(I)+D416*K3(I))
+               CONT(J+NRD*5)=H*(CONT(J+NRD*5)+D513*K4(I)+D514*K10(I)
+     &            +D515*K2(I)+D516*K3(I))
+               CONT(J+NRD*6)=H*(CONT(J+NRD*6)+D613*K4(I)+D614*K10(I)
+     &            +D615*K2(I)+D616*K3(I))
+               CONT(J+NRD*7)=H*(CONT(J+NRD*7)+D713*K4(I)+D714*K10(I)
+     &            +D715*K2(I)+D716*K3(I))
+  63        CONTINUE
+            HOUT=H
+         END IF
+         DO 67 I=1,N
+         K1(I)=K4(I)
+  67     Y(I)=K5(I)
+         XOLD=X
+         X=XPH
+         IF (IOUT.GE.1) THEN
+            CALL SOLOUT(NACCPT+1,XOLD,X,Y,N,CONT,ICOMP,NRD,
+     &                  RPAR,IPAR,IRTRN)
+            IF (IRTRN.LT.0) GOTO 79
+         END IF 
+C ------- NORMAL EXIT
+         IF (LAST) THEN
+            H=HNEW
+            IDID=1
+            RETURN
+         END IF
+         IF(ABS(HNEW).GT.HMAX)HNEW=POSNEG*HMAX  
+         IF(REJECT)HNEW=POSNEG*MIN(ABS(HNEW),ABS(H))
+         REJECT=.FALSE. 
+      ELSE  
+C --- STEP IS REJECTED   
+         HNEW=H/MIN(FACC1,FAC11/SAFE)
+         REJECT=.TRUE.  
+         IF(NACCPT.GE.1)NREJCT=NREJCT+1   
+         LAST=.FALSE.
+      END IF
+      H=HNEW
+      GOTO 1
+C --- FAIL EXIT
+  76  CONTINUE
+      IDID=-4
+      RETURN
+  77  CONTINUE
+      IF (IPRINT.GT.0) WRITE(IPRINT,979)X   
+      IF (IPRINT.GT.0) WRITE(IPRINT,*)' STEP SIZE TOO SMALL, H=',H
+      IDID=-3
+      RETURN
+  78  CONTINUE
+      IF (IPRINT.GT.0) WRITE(IPRINT,979)X   
+      IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &     ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED' 
+      IDID=-2
+      RETURN
+  79  CONTINUE
+      IF (IPRINT.GT.0) WRITE(IPRINT,979)X
+ 979  FORMAT(' EXIT OF DOP853 AT X=',E18.4) 
+      IDID=2
+      RETURN
+      END
+C
+      FUNCTION HINIT853(N,FCN,X,Y,XEND,POSNEG,F0,F1,Y1,IORD,
+     &                       HMAX,ATOL,RTOL,ITOL,RPAR,IPAR)
+C ----------------------------------------------------------
+C ----  COMPUTATION OF AN INITIAL STEP SIZE GUESS
+C ----------------------------------------------------------
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      DIMENSION Y(N),Y1(N),F0(N),F1(N),ATOL(*),RTOL(*)
+      DIMENSION RPAR(*),IPAR(*)
+C ---- COMPUTE A FIRST GUESS FOR EXPLICIT EULER AS
+C ----   H = 0.01 * NORM (Y0) / NORM (F0)
+C ---- THE INCREMENT FOR EXPLICIT EULER IS SMALL
+C ---- COMPARED TO THE SOLUTION
+      DNF=0.0D0
+      DNY=0.0D0 
+      ATOLI=ATOL(1)
+      RTOLI=RTOL(1)    
+      IF (ITOL.EQ.0) THEN   
+        DO 10 I=1,N 
+        SK=ATOLI+RTOLI*ABS(Y(I))
+        DNF=DNF+(F0(I)/SK)**2
+  10    DNY=DNY+(Y(I)/SK)**2 
+      ELSE
+        DO 11 I=1,N 
+        SK=ATOL(I)+RTOL(I)*ABS(Y(I))
+        DNF=DNF+(F0(I)/SK)**2
+  11    DNY=DNY+(Y(I)/SK)**2 
+      END IF
+      IF (DNF.LE.1.D-10.OR.DNY.LE.1.D-10) THEN
+         H=1.0D-6
+      ELSE
+         H=SQRT(DNY/DNF)*0.01D0 
+      END IF
+      H=MIN(H,HMAX)
+      H=SIGN(H,POSNEG) 
+C ---- PERFORM AN EXPLICIT EULER STEP
+      DO 12 I=1,N
+  12  Y1(I)=Y(I)+H*F0(I)
+      CALL FCN(N,X+H,Y1,F1,RPAR,IPAR) 
+C ---- ESTIMATE THE SECOND DERIVATIVE OF THE SOLUTION
+      DER2=0.0D0
+      IF (ITOL.EQ.0) THEN   
+        DO 15 I=1,N 
+        SK=ATOLI+RTOLI*ABS(Y(I))
+  15    DER2=DER2+((F1(I)-F0(I))/SK)**2   
+      ELSE
+        DO 16 I=1,N 
+        SK=ATOL(I)+RTOL(I)*ABS(Y(I))
+  16    DER2=DER2+((F1(I)-F0(I))/SK)**2   
+      END IF
+      DER2=SQRT(DER2)/H
+C ---- STEP SIZE IS COMPUTED SUCH THAT
+C ----  H**IORD * MAX ( NORM (F0), NORM (DER2)) = 0.01
+      DER12=MAX(ABS(DER2),SQRT(DNF))
+      IF (DER12.LE.1.D-15) THEN
+         H1=MAX(1.0D-6,ABS(H)*1.0D-3)
+      ELSE
+         H1=(0.01D0/DER12)**(1.D0/IORD) 
+      END IF
+      H=MIN(100*ABS(H),H1,HMAX)
+      HINIT853=SIGN(H,POSNEG)  
+      RETURN
+      END 
+C
+      FUNCTION CONTD8(II,X,CON,ICOMP,ND)
+C ----------------------------------------------------------
+C     THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT IN CONNECTION
+C     WITH THE OUTPUT-SUBROUTINE FOR DOP853. IT PROVIDES AN
+C     APPROXIMATION TO THE II-TH COMPONENT OF THE SOLUTION AT X.
+C ----------------------------------------------------------
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      DIMENSION CON(8*ND),ICOMP(ND)
+      COMMON /CONDO8/XOLD,H
+C ----- COMPUTE PLACE OF II-TH COMPONENT 
+      I=0 
+      DO 5 J=1,ND 
+      IF (ICOMP(J).EQ.II) I=J
+   5  CONTINUE
+      IF (I.EQ.0) THEN
+         WRITE (6,*) ' NO DENSE OUTPUT AVAILABLE FOR COMP.',II 
+         RETURN
+      END IF  
+      S=(X-XOLD)/H
+      S1=1.D0-S
+      CONPAR=CON(I+ND*4)+S*(CON(I+ND*5)+S1*(CON(I+ND*6)+S*CON(I+ND*7)))
+      CONTD8=CON(I)+S*(CON(I+ND)+S1*(CON(I+ND*2)+S*(CON(I+ND*3)
+     &        +S1*CONPAR)))
+      RETURN
+      END
+

Added: trunk/scipy/integrate/dop/dopri5.f
===================================================================
--- trunk/scipy/integrate/dop/dopri5.f	2009-02-23 05:02:34 UTC (rev 5588)
+++ trunk/scipy/integrate/dop/dopri5.f	2009-02-23 10:25:08 UTC (rev 5589)
@@ -0,0 +1,693 @@
+      SUBROUTINE DOPRI5(N,FCN,X,Y,XEND,
+     &                  RTOL,ATOL,ITOL,
+     &                  SOLOUT,IOUT,
+     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
+C ----------------------------------------------------------
+C     NUMERICAL SOLUTION OF A SYSTEM OF FIRST 0RDER
+C     ORDINARY DIFFERENTIAL EQUATIONS  Y'=F(X,Y).
+C     THIS IS AN EXPLICIT RUNGE-KUTTA METHOD OF ORDER (4)5  
+C     DUE TO DORMAND & PRINCE (WITH STEPSIZE CONTROL AND
+C     DENSE OUTPUT).
+C
+C     AUTHORS: E. HAIRER AND G. WANNER
+C              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
+C              CH-1211 GENEVE 24, SWITZERLAND 
+C              E-MAIL:  Ernst.Hairer@math.unige.ch
+C                       Gerhard.Wanner@math.unige.ch
+C     
+C     THIS CODE IS DESCRIBED IN:
+C         E. HAIRER, S.P. NORSETT AND G. WANNER, SOLVING ORDINARY
+C         DIFFERENTIAL EQUATIONS I. NONSTIFF PROBLEMS. 2ND EDITION.
+C         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
+C         SPRINGER-VERLAG (1993)               
+C      
+C     VERSION OF APRIL 25, 1996
+C     (latest correction of a small bug: August 8, 2005)
+C
+C     INPUT PARAMETERS  
+C     ----------------  
+C     N           DIMENSION OF THE SYSTEM 
+C
+C     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
+C                 VALUE OF F(X,Y):
+C                    SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR)
+C                    DOUBLE PRECISION X,Y(N),F(N)
+C                    F(1)=...   ETC.
+C
+C     X           INITIAL X-VALUE
+C
+C     Y(N)        INITIAL VALUES FOR Y
+C
+C     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
+C
+C     RTOL,ATOL   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
+C                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
+C
+C     ITOL        SWITCH FOR RTOL AND ATOL:
+C                   ITOL=0: BOTH RTOL AND ATOL ARE SCALARS.
+C                     THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF
+C                     Y(I) BELOW RTOL*ABS(Y(I))+ATOL
+C                   ITOL=1: BOTH RTOL AND ATOL ARE VECTORS.
+C                     THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW
+C                     RTOL(I)*ABS(Y(I))+ATOL(I).
+C
+C     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
+C                 NUMERICAL SOLUTION DURING INTEGRATION. 
+C                 IF IOUT.GE.1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
+C                 SUPPLY A DUMMY SUBROUTINE IF IOUT=0. 
+C                 IT MUST HAVE THE FORM
+C                    SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,CON,ICOMP,ND,
+C                                       RPAR,IPAR,IRTRN)
+C                    DIMENSION Y(N),CON(5*ND),ICOMP(ND)
+C                    ....  
+C                 SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
+C                    GRID-POINT "X" (THEREBY THE INITIAL VALUE IS
+C                    THE FIRST GRID-POINT).
+C                 "XOLD" IS THE PRECEEDING GRID-POINT.
+C                 "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
+C                    IS SET <0, DOPRI5 WILL RETURN TO THE CALLING PROGRAM.
+C                    IF THE NUMERICAL SOLUTION IS ALTERED IN SOLOUT,
+C                    SET  IRTRN = 2
+C           
+C          -----  CONTINUOUS OUTPUT: -----
+C                 DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
+C                 FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH
+C                 THE FUNCTION
+C                        >>>   CONTD5(I,S,CON,ICOMP,ND)   <<<
+C                 WHICH PROVIDES AN APPROXIMATION TO THE I-TH
+C                 COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
+C                 S SHOULD LIE IN THE INTERVAL [XOLD,X].
+C
+C     IOUT        SWITCH FOR CALLING THE SUBROUTINE SOLOUT:
+C                    IOUT=0: SUBROUTINE IS NEVER CALLED
+C                    IOUT=1: SUBROUTINE IS USED FOR OUTPUT.
+C                    IOUT=2: DENSE OUTPUT IS PERFORMED IN SOLOUT
+C                            (IN THIS CASE WORK(5) MUST BE SPECIFIED)
+C
+C     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
+C                 WORK(1),...,WORK(20) SERVE AS PARAMETERS FOR THE CODE.
+C                 FOR STANDARD USE, SET THEM TO ZERO BEFORE CALLING.
+C                 "LWORK" MUST BE AT LEAST  8*N+5*NRDENS+21
+C                 WHERE  NRDENS = IWORK(5)
+C
+C     LWORK       DECLARED LENGHT OF ARRAY "WORK".
+C
+C     IWORK       INTEGER WORKING SPACE OF LENGHT "LIWORK".
+C                 IWORK(1),...,IWORK(20) SERVE AS PARAMETERS FOR THE CODE.
+C                 FOR STANDARD USE, SET THEM TO ZERO BEFORE CALLING.
+C                 "LIWORK" MUST BE AT LEAST NRDENS+21 .
+C
+C     LIWORK      DECLARED LENGHT OF ARRAY "IWORK".
+C
+C     RPAR, IPAR  REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH  
+C                 CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING
+C                 PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES. 
+C
+C-----------------------------------------------------------------------
+C 
+C     SOPHISTICATED SETTING OF PARAMETERS
+C     -----------------------------------
+C              SEVERAL PARAMETERS (WORK(1),...,IWORK(1),...) ALLOW
+C              TO ADAPT THE CODE TO THE PROBLEM AND TO THE NEEDS OF
+C              THE USER. FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES.
+C
+C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 2.3D-16.
+C
+C    WORK(2)   THE SAFETY FACTOR IN STEP SIZE PREDICTION,
+C              DEFAULT 0.9D0.
+C
+C    WORK(3), WORK(4)   PARAMETERS FOR STEP SIZE SELECTION
+C              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION
+C                 WORK(3) <= HNEW/HOLD <= WORK(4)
+C              DEFAULT VALUES: WORK(3)=0.2D0, WORK(4)=10.D0
+C
+C    WORK(5)   IS THE "BETA" FOR STABILIZED STEP SIZE CONTROL
+C              (SEE SECTION IV.2). LARGER VALUES OF BETA ( <= 0.1 )
+C              MAKE THE STEP SIZE CONTROL MORE STABLE. DOPRI5 NEEDS
+C              A LARGER BETA THAN HIGHAM & HALL. NEGATIVE WORK(5)
+C              PROVOKE BETA=0.
+C              DEFAULT 0.04D0.
+C
+C    WORK(6)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
+C
+C    WORK(7)   INITIAL STEP SIZE, FOR WORK(7)=0.D0 AN INITIAL GUESS
+C              IS COMPUTED WITH HELP OF THE FUNCTION HINIT
+C
+C    IWORK(1)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
+C              THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000.
+C
+C    IWORK(2)  SWITCH FOR THE CHOICE OF THE COEFFICIENTS
+C              IF IWORK(2).EQ.1  METHOD DOPRI5 OF DORMAND AND PRINCE
+C              (TABLE 5.2 OF SECTION II.5).
+C              AT THE MOMENT THIS IS THE ONLY POSSIBLE CHOICE.
+C              THE DEFAULT VALUE (FOR IWORK(2)=0) IS IWORK(2)=1.
+C
+C    IWORK(3)  SWITCH FOR PRINTING ERROR MESSAGES
+C              IF IWORK(3).LT.0 NO MESSAGES ARE BEING PRINTED
+C              IF IWORK(3).GT.0 MESSAGES ARE PRINTED WITH
+C              WRITE (IWORK(3),*) ...  
+C              DEFAULT VALUE (FOR IWORK(3)=0) IS IWORK(3)=6
+C
+C    IWORK(4)  TEST FOR STIFFNESS IS ACTIVATED AFTER STEP NUMBER
+C              J*IWORK(4) (J INTEGER), PROVIDED IWORK(4).GT.0.
+C              FOR NEGATIVE IWORK(4) THE STIFFNESS TEST IS
+C              NEVER ACTIVATED; DEFAULT VALUE IS IWORK(4)=1000
+C
+C    IWORK(5)  = NRDENS = NUMBER OF COMPONENTS, FOR WHICH DENSE OUTPUT
+C              IS REQUIRED; DEFAULT VALUE IS IWORK(5)=0;
+C              FOR   0 < NRDENS < N   THE COMPONENTS (FOR WHICH DENSE
+C              OUTPUT IS REQUIRED) HAVE TO BE SPECIFIED IN
+C              IWORK(21),...,IWORK(NRDENS+20);
+C              FOR  NRDENS=N  THIS IS DONE BY THE CODE.
+C
+C----------------------------------------------------------------------
+C
+C     OUTPUT PARAMETERS 
+C     ----------------- 
+C     X           X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED
+C                 (AFTER SUCCESSFUL RETURN X=XEND).
+C
+C     Y(N)        NUMERICAL SOLUTION AT X
+C 
+C     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
+C
+C     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
+C                   IDID= 1  COMPUTATION SUCCESSFUL,
+C                   IDID= 2  COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT)
+C                   IDID=-1  INPUT IS NOT CONSISTENT,
+C                   IDID=-2  LARGER NMAX IS NEEDED,
+C                   IDID=-3  STEP SIZE BECOMES TOO SMALL.
+C                   IDID=-4  PROBLEM IS PROBABLY STIFF (INTERRUPTED).
+C
+C   IWORK(17)  NFCN    NUMBER OF FUNCTION EVALUATIONS
+C   IWORK(18)  NSTEP   NUMBER OF COMPUTED STEPS
+C   IWORK(19)  NACCPT  NUMBER OF ACCEPTED STEPS
+C   IWORK(20)  NREJCT  NUMBER OF REJECTED STEPS (DUE TO ERROR TEST),
+C                      (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED)
+C-----------------------------------------------------------------------
+C *** *** *** *** *** *** *** *** *** *** *** *** ***
+C          DECLARATIONS 
+C *** *** *** *** *** *** *** *** *** *** *** *** ***
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      DIMENSION Y(N),ATOL(*),RTOL(*),WORK(LWORK),IWORK(LIWORK)
+      DIMENSION RPAR(*),IPAR(*)
+      LOGICAL ARRET
+      EXTERNAL FCN,SOLOUT
+C *** *** *** *** *** *** ***
+C        SETTING THE PARAMETERS 
+C *** *** *** *** *** *** ***
+      NFCN=0
+      NSTEP=0
+      NACCPT=0
+      NREJCT=0
+      ARRET=.FALSE.
+C -------- IPRINT FOR MONITORING THE PRINTING
+      IF(IWORK(3).EQ.0)THEN
+         IPRINT=6
+      ELSE
+         IPRINT=IWORK(3)
+      END IF
+C -------- NMAX , THE MAXIMAL NUMBER OF STEPS ----- 
+      IF(IWORK(1).EQ.0)THEN
+         NMAX=100000
+      ELSE
+         NMAX=IWORK(1)
+         IF(NMAX.LE.0)THEN
+            IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &          ' WRONG INPUT IWORK(1)=',IWORK(1)
+            ARRET=.TRUE.
+         END IF
+      END IF
+C -------- METH   COEFFICIENTS OF THE METHOD
+      IF(IWORK(2).EQ.0)THEN
+         METH=1
+      ELSE
+         METH=IWORK(2)
+         IF(METH.LE.0.OR.METH.GE.4)THEN
+            IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &          ' CURIOUS INPUT IWORK(2)=',IWORK(2)
+            ARRET=.TRUE.
+         END IF
+      END IF  
+C -------- NSTIFF   PARAMETER FOR STIFFNESS DETECTION  
+      NSTIFF=IWORK(4) 
+      IF (NSTIFF.EQ.0) NSTIFF=1000
+      IF (NSTIFF.LT.0) NSTIFF=NMAX+10
+C -------- NRDENS   NUMBER OF DENSE OUTPUT COMPONENTS
+      NRDENS=IWORK(5)
+      IF(NRDENS.LT.0.OR.NRDENS.GT.N)THEN
+         IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &           ' CURIOUS INPUT IWORK(5)=',IWORK(5)
+         ARRET=.TRUE.
+      ELSE
+            IF(NRDENS.GT.0.AND.IOUT.LT.2)THEN
+               IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &      ' WARNING: PUT IOUT=2 FOR DENSE OUTPUT '
+            END IF 
+            IF (NRDENS.EQ.N) THEN
+                DO 16 I=1,NRDENS
+  16            IWORK(20+I)=I 
+            END IF
+      END IF
+C -------- UROUND   SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0  
+      IF(WORK(1).EQ.0.D0)THEN
+         UROUND=2.3D-16
+      ELSE
+         UROUND=WORK(1)
+         IF(UROUND.LE.1.D-35.OR.UROUND.GE.1.D0)THEN
+            IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &        ' WHICH MACHINE DO YOU HAVE? YOUR UROUND WAS:',WORK(1)
+            ARRET=.TRUE.
+         END IF
+      END IF
+C -------  SAFETY FACTOR -------------
+      IF(WORK(2).EQ.0.D0)THEN
+         SAFE=0.9D0
+      ELSE
+         SAFE=WORK(2)
+         IF(SAFE.GE.1.D0.OR.SAFE.LE.1.D-4)THEN
+            IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &          ' CURIOUS INPUT FOR SAFETY FACTOR WORK(2)=',WORK(2)
+            ARRET=.TRUE.
+         END IF
+      END IF
+C -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
+      IF(WORK(3).EQ.0.D0)THEN
+         FAC1=0.2D0
+      ELSE
+         FAC1=WORK(3)
+      END IF
+      IF(WORK(4).EQ.0.D0)THEN
+         FAC2=10.D0
+      ELSE
+         FAC2=WORK(4)
+      END IF
+C --------- BETA FOR STEP CONTROL STABILIZATION -----------
+      IF(WORK(5).EQ.0.D0)THEN
+         BETA=0.04D0
+      ELSE
+         IF(WORK(5).LT.0.D0)THEN
+            BETA=0.D0
+         ELSE
+            BETA=WORK(5)
+            IF(BETA.GT.0.2D0)THEN
+               IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &          ' CURIOUS INPUT FOR BETA: WORK(5)=',WORK(5)
+            ARRET=.TRUE.
+         END IF
+         END IF
+      END IF
+C -------- MAXIMAL STEP SIZE
+      IF(WORK(6).EQ.0.D0)THEN
+         HMAX=XEND-X
+      ELSE
+         HMAX=WORK(6)
+      END IF
+C -------- INITIAL STEP SIZE
+      H=WORK(7)
+C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
+      IEY1=21
+      IEK1=IEY1+N
+      IEK2=IEK1+N
+      IEK3=IEK2+N
+      IEK4=IEK3+N
+      IEK5=IEK4+N
+      IEK6=IEK5+N
+      IEYS=IEK6+N
+      IECO=IEYS+N
+C ------ TOTAL STORAGE REQUIREMENT -----------
+      ISTORE=IEYS+5*NRDENS-1
+      IF(ISTORE.GT.LWORK)THEN
+        IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &   ' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
+        ARRET=.TRUE.
+      END IF
+      ICOMP=21
+      ISTORE=ICOMP+NRDENS-1
+      IF(ISTORE.GT.LIWORK)THEN
+        IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &   ' INSUFFICIENT STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
+        ARRET=.TRUE.
+      END IF
+C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
+      IF (ARRET) THEN
+         IDID=-1
+         RETURN
+      END IF
+C -------- CALL TO CORE INTEGRATOR ------------
+      CALL DOPCOR(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,IPRINT,
+     &   SOLOUT,IOUT,IDID,NMAX,UROUND,METH,NSTIFF,SAFE,BETA,FAC1,FAC2,
+     &   WORK(IEY1),WORK(IEK1),WORK(IEK2),WORK(IEK3),WORK(IEK4),
+     &   WORK(IEK5),WORK(IEK6),WORK(IEYS),WORK(IECO),IWORK(ICOMP),
+     &   NRDENS,RPAR,IPAR,NFCN,NSTEP,NACCPT,NREJCT)
+      WORK(7)=H
+      IWORK(17)=NFCN
+      IWORK(18)=NSTEP
+      IWORK(19)=NACCPT
+      IWORK(20)=NREJCT
+C ----------- RETURN -----------
+      RETURN
+      END
+C
+C
+C
+C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
+C
+      SUBROUTINE DOPCOR(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,IPRINT,
+     &   SOLOUT,IOUT,IDID,NMAX,UROUND,METH,NSTIFF,SAFE,BETA,FAC1,FAC2,
+     &   Y1,K1,K2,K3,K4,K5,K6,YSTI,CONT,ICOMP,NRD,RPAR,IPAR,
+     &   NFCN,NSTEP,NACCPT,NREJCT)
+C ----------------------------------------------------------
+C     CORE INTEGRATOR FOR DOPRI5
+C     PARAMETERS SAME AS IN DOPRI5 WITH WORKSPACE ADDED 
+C ---------------------------------------------------------- 
+C         DECLARATIONS 
+C ---------------------------------------------------------- 
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      DOUBLE PRECISION K1(N),K2(N),K3(N),K4(N),K5(N),K6(N)
+      DIMENSION Y(N),Y1(N),YSTI(N),ATOL(*),RTOL(*),RPAR(*),IPAR(*)
+      DIMENSION CONT(5*NRD),ICOMP(NRD)
+      LOGICAL REJECT,LAST 
+      EXTERNAL FCN
+      COMMON /CONDO5/XOLD,HOUT
+C *** *** *** *** *** *** ***
+C  INITIALISATIONS
+C *** *** *** *** *** *** *** 
+      IF (METH.EQ.1) CALL CDOPRI(C2,C3,C4,C5,E1,E3,E4,E5,E6,E7,
+     &                    A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,
+     &                    A61,A62,A63,A64,A65,A71,A73,A74,A75,A76,
+     &                    D1,D3,D4,D5,D6,D7)
+      FACOLD=1.D-4  
+      EXPO1=0.2D0-BETA*0.75D0
+      FACC1=1.D0/FAC1
+      FACC2=1.D0/FAC2
+      POSNEG=SIGN(1.D0,XEND-X)  
+C --- INITIAL PREPARATIONS   
+      ATOLI=ATOL(1)
+      RTOLI=RTOL(1)    
+      LAST=.FALSE. 
+      HLAMB=0.D0
+      IASTI=0
+      CALL FCN(N,X,Y,K1,RPAR,IPAR)
+      HMAX=ABS(HMAX)     
+      IORD=5  
+      IF (H.EQ.0.D0) H=HINIT(N,FCN,X,Y,XEND,POSNEG,K1,K2,K3,IORD,
+     &                       HMAX,ATOL,RTOL,ITOL,RPAR,IPAR)
+      NFCN=NFCN+2
+      REJECT=.FALSE.
+      XOLD=X
+      IF (IOUT.NE.0) THEN 
+          IRTRN=1
+          HOUT=H
+          CALL SOLOUT(NACCPT+1,XOLD,X,Y,N,CONT,ICOMP,NRD,
+     &                RPAR,IPAR,IRTRN)
+          IF (IRTRN.LT.0) GOTO 79
+      ELSE
+          IRTRN=0
+      END IF
+C --- BASIC INTEGRATION STEP  
+   1  CONTINUE
+      IF (NSTEP.GT.NMAX) GOTO 78
+      IF (0.1D0*ABS(H).LE.ABS(X)*UROUND)GOTO 77
+      IF ((X+1.01D0*H-XEND)*POSNEG.GT.0.D0) THEN
+         H=XEND-X 
+         LAST=.TRUE.
+      END IF
+      NSTEP=NSTEP+1
+C --- THE FIRST 6 STAGES
+      IF (IRTRN.GE.2) THEN
+         CALL FCN(N,X,Y,K1,RPAR,IPAR)
+      END IF
+      DO 22 I=1,N 
+  22  Y1(I)=Y(I)+H*A21*K1(I)
+      CALL FCN(N,X+C2*H,Y1,K2,RPAR,IPAR)
+      DO 23 I=1,N 
+  23  Y1(I)=Y(I)+H*(A31*K1(I)+A32*K2(I))
+      CALL FCN(N,X+C3*H,Y1,K3,RPAR,IPAR)
+      DO 24 I=1,N 
+  24  Y1(I)=Y(I)+H*(A41*K1(I)+A42*K2(I)+A43*K3(I))
+      CALL FCN(N,X+C4*H,Y1,K4,RPAR,IPAR)
+      DO 25 I=1,N 
+  25  Y1(I)=Y(I)+H*(A51*K1(I)+A52*K2(I)+A53*K3(I)+A54*K4(I))
+      CALL FCN(N,X+C5*H,Y1,K5,RPAR,IPAR) 
+      DO 26 I=1,N 
+  26  YSTI(I)=Y(I)+H*(A61*K1(I)+A62*K2(I)+A63*K3(I)+A64*K4(I)+A65*K5(I))
+      XPH=X+H
+      CALL FCN(N,XPH,YSTI,K6,RPAR,IPAR)
+      DO 27 I=1,N 
+  27  Y1(I)=Y(I)+H*(A71*K1(I)+A73*K3(I)+A74*K4(I)+A75*K5(I)+A76*K6(I))  
+      CALL FCN(N,XPH,Y1,K2,RPAR,IPAR)
+      IF (IOUT.GE.2) THEN 
+            DO 40 J=1,NRD
+            I=ICOMP(J)
+            CONT(4*NRD+J)=H*(D1*K1(I)+D3*K3(I)+D4*K4(I)+D5*K5(I)
+     &                   +D6*K6(I)+D7*K2(I)) 
+  40        CONTINUE
+      END IF
+      DO 28 I=1,N 
+  28  K4(I)=(E1*K1(I)+E3*K3(I)+E4*K4(I)+E5*K5(I)+E6*K6(I)+E7*K2(I))*H
+      NFCN=NFCN+6 
+C --- ERROR ESTIMATION  
+      ERR=0.D0
+      IF (ITOL.EQ.0) THEN   
+        DO 41 I=1,N 
+        SK=ATOLI+RTOLI*MAX(ABS(Y(I)),ABS(Y1(I)))
+  41    ERR=ERR+(K4(I)/SK)**2
+      ELSE
+        DO 42 I=1,N 
+        SK=ATOL(I)+RTOL(I)*MAX(ABS(Y(I)),ABS(Y1(I)))
+  42    ERR=ERR+(K4(I)/SK)**2
+      END IF
+      ERR=SQRT(ERR/N)  
+C --- COMPUTATION OF HNEW
+      FAC11=ERR**EXPO1
+C --- LUND-STABILIZATION
+      FAC=FAC11/FACOLD**BETA
+C --- WE REQUIRE  FAC1 <= HNEW/H <= FAC2
+      FAC=MAX(FACC2,MIN(FACC1,FAC/SAFE))
+      HNEW=H/FAC  
+      IF(ERR.LE.1.D0)THEN
+C --- STEP IS ACCEPTED  
+         FACOLD=MAX(ERR,1.0D-4)
+         NACCPT=NACCPT+1
+C ------- STIFFNESS DETECTION
+         IF (MOD(NACCPT,NSTIFF).EQ.0.OR.IASTI.GT.0) THEN
+            STNUM=0.D0
+            STDEN=0.D0
+            DO 64 I=1,N 
+               STNUM=STNUM+(K2(I)-K6(I))**2
+               STDEN=STDEN+(Y1(I)-YSTI(I))**2
+ 64         CONTINUE  
+            IF (STDEN.GT.0.D0) HLAMB=H*SQRT(STNUM/STDEN) 
+            IF (HLAMB.GT.3.25D0) THEN
+               NONSTI=0
+               IASTI=IASTI+1  
+               IF (IASTI.EQ.15) THEN
+                  IF (IPRINT.GT.0) WRITE (IPRINT,*) 
+     &               ' THE PROBLEM SEEMS TO BECOME STIFF AT X = ',X   
+                  IF (IPRINT.LE.0) GOTO 76
+               END IF
+            ELSE
+               NONSTI=NONSTI+1  
+               IF (NONSTI.EQ.6) IASTI=0
+            END IF
+         END IF 
+         IF (IOUT.GE.2) THEN 
+            DO 43 J=1,NRD
+            I=ICOMP(J)
+            YD0=Y(I)
+            YDIFF=Y1(I)-YD0
+            BSPL=H*K1(I)-YDIFF 
+            CONT(J)=Y(I)
+            CONT(NRD+J)=YDIFF
+            CONT(2*NRD+J)=BSPL
+            CONT(3*NRD+J)=-H*K2(I)+YDIFF-BSPL
+  43        CONTINUE
+         END IF
+         DO 44 I=1,N
+         K1(I)=K2(I)
+  44     Y(I)=Y1(I)
+         XOLD=X
+         X=XPH
+         IF (IOUT.NE.0) THEN
+            HOUT=H
+            CALL SOLOUT(NACCPT+1,XOLD,X,Y,N,CONT,ICOMP,NRD,
+     &                  RPAR,IPAR,IRTRN)
+            IF (IRTRN.LT.0) GOTO 79
+         END IF 
+C ------- NORMAL EXIT
+         IF (LAST) THEN
+            H=HNEW
+            IDID=1
+            RETURN
+         END IF
+         IF(ABS(HNEW).GT.HMAX)HNEW=POSNEG*HMAX  
+         IF(REJECT)HNEW=POSNEG*MIN(ABS(HNEW),ABS(H))
+         REJECT=.FALSE. 
+      ELSE  
+C --- STEP IS REJECTED   
+         HNEW=H/MIN(FACC1,FAC11/SAFE)
+         REJECT=.TRUE.  
+         IF(NACCPT.GE.1)NREJCT=NREJCT+1   
+         LAST=.FALSE.
+      END IF
+      H=HNEW
+      GOTO 1
+C --- FAIL EXIT
+  76  CONTINUE
+      IDID=-4
+      RETURN
+  77  CONTINUE
+      IF (IPRINT.GT.0) WRITE(IPRINT,979)X   
+      IF (IPRINT.GT.0) WRITE(IPRINT,*)' STEP SIZE T0O SMALL, H=',H
+      IDID=-3
+      RETURN
+  78  CONTINUE
+      IF (IPRINT.GT.0) WRITE(IPRINT,979)X   
+      IF (IPRINT.GT.0) WRITE(IPRINT,*)
+     &     ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED' 
+      IDID=-2
+      RETURN
+  79  CONTINUE
+      IF (IPRINT.GT.0) WRITE(IPRINT,979)X
+ 979  FORMAT(' EXIT OF DOPRI5 AT X=',E18.4) 
+      IDID=2
+      RETURN
+      END
+C
+      FUNCTION HINIT(N,FCN,X,Y,XEND,POSNEG,F0,F1,Y1,IORD,
+     &                 HMAX,ATOL,RTOL,ITOL,RPAR,IPAR)
+C ----------------------------------------------------------
+C ----  COMPUTATION OF AN INITIAL STEP SIZE GUESS
+C ----------------------------------------------------------
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      DIMENSION Y(N),Y1(N),F0(N),F1(N),ATOL(*),RTOL(*)
+      DIMENSION RPAR(*),IPAR(*)
+C ---- COMPUTE A FIRST GUESS FOR EXPLICIT EULER AS
+C ----   H = 0.01 * NORM (Y0) / NORM (F0)
+C ---- THE INCREMENT FOR EXPLICIT EULER IS SMALL
+C ---- COMPARED TO THE SOLUTION
+      DNF=0.0D0
+      DNY=0.0D0 
+      ATOLI=ATOL(1)
+      RTOLI=RTOL(1)    
+      IF (ITOL.EQ.0) THEN   
+        DO 10 I=1,N 
+        SK=ATOLI+RTOLI*ABS(Y(I))
+        DNF=DNF+(F0(I)/SK)**2
+  10    DNY=DNY+(Y(I)/SK)**2 
+      ELSE
+        DO 11 I=1,N 
+        SK=ATOL(I)+RTOL(I)*ABS(Y(I))
+        DNF=DNF+(F0(I)/SK)**2
+  11    DNY=DNY+(Y(I)/SK)**2 
+      END IF
+      IF (DNF.LE.1.D-10.OR.DNY.LE.1.D-10) THEN
+         H=1.0D-6
+      ELSE
+         H=SQRT(DNY/DNF)*0.01D0 
+      END IF
+      H=MIN(H,HMAX)
+      H=SIGN(H,POSNEG) 
+C ---- PERFORM AN EXPLICIT EULER STEP
+      DO 12 I=1,N
+  12  Y1(I)=Y(I)+H*F0(I)
+      CALL FCN(N,X+H,Y1,F1,RPAR,IPAR) 
+C ---- ESTIMATE THE SECOND DERIVATIVE OF THE SOLUTION
+      DER2=0.0D0
+      IF (ITOL.EQ.0) THEN   
+        DO 15 I=1,N 
+        SK=ATOLI+RTOLI*ABS(Y(I))
+  15    DER2=DER2+((F1(I)-F0(I))/SK)**2   
+      ELSE
+        DO 16 I=1,N 
+        SK=ATOL(I)+RTOL(I)*ABS(Y(I))
+  16    DER2=DER2+((F1(I)-F0(I))/SK)**2   
+      END IF
+      DER2=SQRT(DER2)/H
+C ---- STEP SIZE IS COMPUTED SUCH THAT
+C ----  H**IORD * MAX ( NORM (F0), NORM (DER2)) = 0.01
+      DER12=MAX(ABS(DER2),SQRT(DNF))
+      IF (DER12.LE.1.D-15) THEN
+         H1=MAX(1.0D-6,ABS(H)*1.0D-3)
+      ELSE
+         H1=(0.01D0/DER12)**(1.D0/IORD) 
+      END IF
+      H=MIN(100*ABS(H),H1,HMAX)
+      HINIT=SIGN(H,POSNEG)  
+      RETURN
+      END 
+C
+      FUNCTION CONTD5(II,X,CON,ICOMP,ND)
+C ----------------------------------------------------------
+C     THIS FUNCTION CAN BE USED FOR CONTINUOUS OUTPUT IN CONNECTION
+C     WITH THE OUTPUT-SUBROUTINE FOR DOPRI5. IT PROVIDES AN
+C     APPROXIMATION TO THE II-TH COMPONENT OF THE SOLUTION AT X.
+C ----------------------------------------------------------
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      DIMENSION CON(5*ND),ICOMP(ND)
+      COMMON /CONDO5/XOLD,H
+C ----- COMPUTE PLACE OF II-TH COMPONENT 
+      I=0 
+      DO 5 J=1,ND 
+      IF (ICOMP(J).EQ.II) I=J
+   5  CONTINUE
+      IF (I.EQ.0) THEN
+         WRITE (6,*) ' NO DENSE OUTPUT AVAILABLE FOR COMP.',II 
+         RETURN
+      END IF  
+      THETA=(X-XOLD)/H
+      THETA1=1.D0-THETA
+      CONTD5=CON(I)+THETA*(CON(ND+I)+THETA1*(CON(2*ND+I)+THETA*
+     &           (CON(3*ND+I)+THETA1*CON(4*ND+I))))
+      RETURN
+      END
+C
+      SUBROUTINE CDOPRI(C2,C3,C4,C5,E1,E3,E4,E5,E6,E7,
+     &                    A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,
+     &                    A61,A62,A63,A64,A65,A71,A73,A74,A75,A76,
+     &                    D1,D3,D4,D5,D6,D7)
+C ----------------------------------------------------------
+C     RUNGE-KUTTA COEFFICIENTS OF DORMAND AND PRINCE (1980)
+C ----------------------------------------------------------
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      C2=0.2D0
+      C3=0.3D0
+      C4=0.8D0
+      C5=8.D0/9.D0
+      A21=0.2D0
+      A31=3.D0/40.D0
+      A32=9.D0/40.D0
+      A41=44.D0/45.D0
+      A42=-56.D0/15.D0
+      A43=32.D0/9.D0
+      A51=19372.D0/6561.D0
+      A52=-25360.D0/2187.D0
+      A53=64448.D0/6561.D0
+      A54=-212.D0/729.D0
+      A61=9017.D0/3168.D0
+      A62=-355.D0/33.D0
+      A63=46732.D0/5247.D0
+      A64=49.D0/176.D0
+      A65=-5103.D0/18656.D0
+      A71=35.D0/384.D0
+      A73=500.D0/1113.D0
+      A74=125.D0/192.D0
+      A75=-2187.D0/6784.D0
+      A76=11.D0/84.D0
+      E1=71.D0/57600.D0
+      E3=-71.D0/16695.D0
+      E4=71.D0/1920.D0
+      E5=-17253.D0/339200.D0
+      E6=22.D0/525.D0
+      E7=-1.D0/40.D0  
+C ---- DENSE OUTPUT OF SHAMPINE (1986)
+      D1=-12715105075.D0/11282082432.D0
+      D3=87487479700.D0/32700410799.D0
+      D4=-10690763975.D0/1880347072.D0
+      D5=701980252875.D0/199316789632.D0
+      D6=-1453857185.D0/822651844.D0
+      D7=69997945.D0/29380423.D0
+      RETURN
+      END
+

Added: trunk/scipy/integrate/dop.pyf
===================================================================
--- trunk/scipy/integrate/dop.pyf	2009-02-23 05:02:34 UTC (rev 5588)
+++ trunk/scipy/integrate/dop.pyf	2009-02-23 10:25:08 UTC (rev 5589)
@@ -0,0 +1,80 @@
+!%f90 -*- f90 -*-
+!Author: John Travers
+!Date:   22 Feb 2009
+
+python module __user__routines 
+    interface
+       subroutine fcn(n,x,y,f,rpar,ipar)
+         integer intent(hide) :: n
+         double precision intent(in) :: x
+         double precision dimension(n),intent(in,c) :: y
+         double precision dimension(n),intent(out,c) :: f
+         double precision intent(hide) :: rpar
+         integer intent(hide) :: ipar
+       end subroutine fcn
+       subroutine solout(nr,xold,x,y,n,con,icomp,nd,rpar,ipar,irtn)
+         integer intent(in) :: nr
+         integer intent(hide) :: n
+         double precision intent(in) :: xold, x
+         double precision dimension(n),intent(c,in) :: y
+         integer intent(in) :: nd
+         integer dimension(nd), intent(in) :: icomp
+         double precision dimension(5*nd), intent(in) :: con
+         double precision intent(hide) :: rpar
+         integer intent(hide) :: ipar
+         integer intent(out) :: irtn
+       end subroutine solout
+    end interface
+end python module __user__routines
+
+python module dop
+    interface
+       subroutine dopri5(n,fcn,x,y,xend,rtol,atol,itol,solout,iout,work,lwork,iwork,liwork,rpar,ipar,idid)
+         use __user__routines
+         external fcn
+         external solout
+         integer intent(hide),depend(y) :: n = len(y)
+         double precision dimension(n),intent(in,out,copy) :: y
+         double precision intent(in,out):: x
+         double precision intent(in):: xend
+         double precision dimension(*),intent(in),check(len(atol)<&
+              &=1||len(atol)>=n),depend(n) :: atol
+         double precision dimension(*),intent(in),check(len(rtol)==len(atol)), &
+              depend(atol) :: rtol
+         integer intent(hide), depend(atol) :: itol = (len(atol)<=1?0:1)
+         integer intent(hide) :: iout=0
+         double precision dimension(*), intent(in), check(len(work)>=8*n+21), &
+                              :: work
+         integer intent(hide), depend(work) :: lwork = len(work)
+         integer intent(in,out), dimension(*), check(len(iwork)>=21) :: iwork
+         integer intent(hide), depend(iwork) :: liwork = len(iwork)
+         integer intent(out) :: idid
+         double precision intent(hide) :: rpar = 0.0
+         integer intent(hide) :: ipar = 0
+       end subroutine dopri5
+       subroutine dop853(n,fcn,x,y,xend,rtol,atol,itol,solout,iout,work,lwork,iwork,liwork,rpar,ipar,idid)
+         use __user__routines
+         external fcn
+         external solout
+         integer intent(hide),depend(y) :: n = len(y)
+         double precision dimension(n),intent(in,out,copy) :: y
+         double precision intent(in,out):: x
+         double precision intent(in):: xend
+         double precision dimension(*),intent(in),check(len(atol)<&
+              &=1||len(atol)>=n),depend(n) :: atol
+         double precision dimension(*),intent(in),check(len(rtol)==len(atol)), &
+              depend(atol) :: rtol
+         integer intent(hide), depend(atol) :: itol = (len(atol)<=1?0:1)
+         integer intent(hide) :: iout=0
+         double precision dimension(*), intent(in), check(len(work)>=8*n+21), &
+                              :: work
+         integer intent(hide), depend(work) :: lwork = len(work)
+         integer intent(in,out), dimension(*), check(len(iwork)>=21) :: iwork
+         integer intent(hide), depend(iwork) :: liwork = len(iwork)
+         integer intent(out) :: idid
+         double precision intent(hide) :: rpar = 0.0
+         integer intent(hide) :: ipar = 0
+       end subroutine dop853
+    end interface
+end python module dop
+

Modified: trunk/scipy/integrate/ode.py
===================================================================
--- trunk/scipy/integrate/ode.py	2009-02-23 05:02:34 UTC (rev 5588)
+++ trunk/scipy/integrate/ode.py	2009-02-23 10:25:08 UTC (rev 5589)
@@ -1,4 +1,4 @@
-# Authors: Pearu Peterson, Pauli Virtanen
+# Authors: Pearu Peterson, Pauli Virtanen, John Travers
 """
 First-order ODE integrators
 
@@ -97,6 +97,58 @@
     failures, and for this problem one should instead use DVODE on the
     equivalent real system (in the real and imaginary parts of y).
 
+dopri5
+~~~~~~
+
+    Numerical solution of a system of first order
+    ordinary differential equations  y'=f(x,y).
+    this is an explicit runge-kutta method of order (4)5  
+    due to Dormand & Prince (with stepsize control and
+    dense output).
+
+    Authors: E. Hairer and G. Wanner
+             Universite de Geneve, Dept. de Mathematiques
+             CH-1211 Geneve 24, Switzerland 
+             e-mail:  ernst.hairer@math.unige.ch
+                      gerhard.wanner@math.unige.ch
+     
+    This code is described in: 
+          E. Hairer, S.P. Norsett and G. Wanner, Solving Ordinary
+          Differential Equations i. Nonstiff Problems. 2nd edition.
+          Springer Series in Computational Mathematics,
+          Springer-Verlag (1993)               
+
+This integrator accepts the following parameters in set_integrator()
+method of the ode class:
+
+- atol : float or sequence
+  absolute tolerance for solution
+- rtol : float or sequence
+  relative tolerance for solution
+- nsteps : int
+  Maximum number of (internally defined) steps allowed during one
+  call to the solver.
+- first_step : float
+- max_step : float
+- safety : float
+  Safety factor on new step selection (default 0.9)
+- ifactor : float
+- dfactor : float
+  Maximum factor to increase/decrease step sixe by in one step
+- beta : float
+  Beta parameter for stabilised step size control.
+
+dop853
+~~~~~~
+
+    Numerical solution of a system of first 0rder
+    ordinary differential equations  y'=f(x,y).
+    this is an explicit runge-kutta method of order 8(5,3)  
+    due to Dormand & Prince (with stepsize control and
+    dense output).
+
+    Options and references the same as dopri5.
+
 """
 
 if __doc__:
@@ -153,6 +205,7 @@
 from numpy import asarray, array, zeros, int32, isscalar
 
 import vode as _vode
+import dop as _dop
 
 #------------------------------------------------------------------------------
 # User interface
@@ -564,3 +617,110 @@
 
 if zvode.runner is not None:
     IntegratorBase.integrator_classes.append(zvode)
+
+class dopri5(IntegratorBase):
+
+    runner = getattr(_dop,'dopri5',None)
+
+    messages = { 1 : 'computation successful',
+                 2 : 'comput. successful (interrupted by solout)',
+                -1 : 'input is not consistent',
+                -2 : 'larger nmax is needed',
+                -3 : 'step size becomes too small',
+                -4 : 'problem is probably stiff (interrupted)',
+               }
+
+    def __init__(self,
+                 rtol=1e-6,atol=1e-12,
+                 nsteps = 500,
+                 max_step = 0.0,
+                 first_step = 0.0, # determined by solver
+                 safety = 0.9,
+                 ifactor = 10.0,
+                 dfactor = 0.2,
+                 beta = 0.0,
+                 method = None
+                 ):
+        self.rtol = rtol
+        self.atol = atol
+        self.nsteps = nsteps
+        self.max_step = max_step
+        self.first_step = first_step
+        self.safety = safety
+        self.ifactor = ifactor
+        self.dfactor = dfactor
+        self.beta = beta
+        self.success = 1
+
+    def reset(self,n,has_jac):
+        work = zeros((8*n+21,), float)
+        work[1] = self.safety
+        work[2] = self.dfactor
+        work[3] = self.ifactor
+        work[4] = self.beta
+        work[5] = self.max_step
+        work[6] = self.first_step
+        self.work = work
+        iwork = zeros((21,), int32)
+        iwork[0] = self.nsteps
+        self.iwork = iwork
+        self.call_args = [self.rtol,self.atol,self._solout,self.work,self.iwork]
+        self.success = 1
+
+    def run(self,f,jac,y0,t0,t1,f_params,jac_params):
+        x,y,iwork,idid = self.runner(*((f,t0,y0,t1) + tuple(self.call_args)))
+        if idid < 0:
+            print 'dopri5:',self.messages.get(idid,'Unexpected idid=%s'%idid)
+            self.success = 0
+        return y,x
+        
+    def _solout(self, *args):
+        # dummy solout function
+        pass
+
+if dopri5.runner:
+    IntegratorBase.integrator_classes.append(dopri5)
+
+class dop853(dopri5):
+
+    runner = getattr(_dop,'dop853',None)
+
+    def __init__(self,
+                 rtol=1e-6,atol=1e-12,
+                 nsteps = 500,
+                 max_step = 0.0,
+                 first_step = 0.0, # determined by solver
+                 safety = 0.9,
+                 ifactor = 6.0,
+                 dfactor = 0.3,
+                 beta = 0.0,
+                 method = None
+                 ):
+        self.rtol = rtol
+        self.atol = atol
+        self.nsteps = nsteps
+        self.max_step = max_step
+        self.first_step = first_step
+        self.safety = safety
+        self.ifactor = ifactor
+        self.dfactor = dfactor
+        self.beta = beta
+        self.success = 1
+
+    def reset(self,n,has_jac):
+        work = zeros((11*n+21,), float)
+        work[1] = self.safety
+        work[2] = self.dfactor
+        work[3] = self.ifactor
+        work[4] = self.beta
+        work[5] = self.max_step
+        work[6] = self.first_step
+        self.work = work
+        iwork = zeros((21,), int32)
+        iwork[0] = self.nsteps
+        self.iwork = iwork
+        self.call_args = [self.rtol,self.atol,self._solout,self.work,self.iwork]
+        self.success = 1
+
+if dop853.runner:
+    IntegratorBase.integrator_classes.append(dop853)

Modified: trunk/scipy/integrate/setup.py
===================================================================
--- trunk/scipy/integrate/setup.py	2009-02-23 05:02:34 UTC (rev 5588)
+++ trunk/scipy/integrate/setup.py	2009-02-23 10:25:08 UTC (rev 5589)
@@ -18,6 +18,8 @@
                        sources=[join('quadpack','*.f')])
     config.add_library('odepack',
                        sources=[join('odepack','*.f')])
+    config.add_library('dop',
+                       sources=[join('dop','*.f')])
     # should we try to weed through files and replace with calls to
     # LAPACK routines?
     # Yes, someday...
@@ -33,6 +35,7 @@
                          depends=['quadpack.h','__quadpack.h'])
     # odepack
     libs = ['odepack','linpack_lite','mach']
+    
 
     # Remove libraries key from blas_opt
     if 'libraries' in blas_opt:    # key doesn't exist on OS X ...
@@ -54,6 +57,11 @@
                          libraries=libs,
                          **newblas)
 
+    # dop
+    config.add_extension('dop',
+                         sources=['dop.pyf'],
+                         libraries=['dop'])
+
     config.add_data_dir('tests')
     return config
 

Modified: trunk/scipy/integrate/tests/test_integrate.py
===================================================================
--- trunk/scipy/integrate/tests/test_integrate.py	2009-02-23 05:02:34 UTC (rev 5588)
+++ trunk/scipy/integrate/tests/test_integrate.py	2009-02-23 10:25:08 UTC (rev 5589)
@@ -69,6 +69,22 @@
                 self._do_problem(problem, 'zvode', 'adams')
             self._do_problem(problem, 'zvode', 'bdf')
 
+    def test_dopri5(self):
+        """Check the dopri5 solver"""
+        for problem_cls in PROBLEMS:
+            problem = problem_cls()
+            if problem.cmplx: continue
+            if problem.stiff: continue
+            self._do_problem(problem, 'dopri5')
+            
+    def test_dop853(self):
+        """Check the dop853 solver"""
+        for problem_cls in PROBLEMS:
+            problem = problem_cls()
+            if problem.cmplx: continue
+            if problem.stiff: continue
+            self._do_problem(problem, 'dop853')
+
 #------------------------------------------------------------------------------
 # Test problems
 #------------------------------------------------------------------------------



More information about the Scipy-svn mailing list