C * * * * * * * * * * * * * * * * * * * * * * * * *
C --- DRIVER FOR MEBDF ON THE ROBERTSON PROBLEM
C * * * * * * * * * * * * * * * * * * * * * * * * *
        IMPLICIT REAL*8 (A-H,O-Z)
C --- PARAMETERS FOR LSODE (FULL JACOBIAN)
        PARAMETER (ND=3,LWORK=(32+3*ND)*ND+3,LIWORK=ND+14)
        DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd),
     +    error(nd),massbnd(4)
        REAL*4 TARRAY(2)
        External f,pderv,mas
c ------ FILE DE DONNEES ----------
        write(6,2050)
2050  format(1x,'mebdf results for robertson')
C --- LOOP FOR DIFFERENT TOLERANCES
        NTOLMN=2
        NTOLMX=10
        NTOLDF=4
        NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1
        TOLST=0.1D0**NTOLMN
        TOLFC=0.1D0**(1.D0/NTOLDF)
        iwork(14)=100000
        work(1)=0.0D+0
        DO 30 NTOL=1,NRLOOP
        summ=0.0D+0
      sumh=0.0D+0
C --- DIMENSION OF THE SYSTEM
        N=3
C --- SET DEFAULT VALUES 
        MF=21 
        INDEX=1
        LOUT=6
C --- INITIAL VALUES
        X=0.0D0
        Y(1)=1.0D0
        Y(2)=0.0D0
        y(3)=0.0D+0
C --- REQUIRED TOLERANCE
        RTOL=TOLST
        ATOL=RTOL*1.0D-6
        ITOL=2
        MAXDER=6
C --- MAXIMAL NUMBER OF STEPS
c        IWORK(6)=10000
C --- ENDPOINT OF INTEGRATION
        XEND=1.0d+11
        H=1.0D-6
         it1=mclock()
        XOUT=1.0D+0
        massbnd(1)=0
        DO 20 I=1,12
 220    CONTINUE
C --- CALL OF THE SUBROUTINE  
        CALL mebdf(N,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, 
     +       WORK,LIWORK,IWORK,MBND,massbnd,MAXDER,ITOL,RTOL,ATOL,
     +       RPAR,IPAR,F,PDERV,MAS,ierr)
        IF (INDEX.EQ.1) THEN
           INDEX=0
           GOTO 220
        ELSE
c         if (mod(ntol,4).eq.1) write (6,*) xout,y(1),y(2)
         do 151 k=1,3
c         write(6,*) y(k)
151      continue
c          write(6,*) h
           XOUT=XOUT*10.0D+0
        END IF

      if(i.eq.1) then
      true(1)= 0.9664597373330035E+00
       true(2)=0.3074626578578675E-04
      true(3)=0.3350951640121071E-01
      else if(i.eq.2) then
      true(1)=0.8413699238414729E+00
      true(2)=0.1623390937990473E-04
      true(3)=0.1586138422491472E+00
      else if(i.eq.3) then
      true(1)=0.6172348823960878E+00
       true(2)=0.6153591274639123E-05
      true(3)=0.3827589640126376E+00
      else if(i.eq.4) then
      true(1)=0.3368745306607069E+00
      true(2)=0.2013702318261393E-05
      true(3)=0.6631234556369748E+00
      else if(i.eq.5) then
      true(1)=0.1073004285378040E+00
      true(2)=0.4800166972571660E-06
      true(3)=0.8926990914454987E+00
      else if(i.eq.6) then
      true(1)=0.1786592114209946E-01
      true(2)=0.7274751468436319E-07
      true(3)=0.9821340061103859E+00
      else if(i.eq.7) then
      true(1)= 0.2031483924973415E-02
       true(2)=0.8142277783356159E-08
      true(3)=0.9979685079327488E+00
      else if(i.eq.8) then
      true(1)=0.2076093439016395E-03
      true(2)=0.8306077485067610E-09
      true(3)=0.9997923898254906E+00
      else if(i.eq.9) then
      true(1)=0.2082417512179460E-04
      true(2)=0.8329841429908955E-10
      true(3)=0.9999791757415798E+00
      else if(i.eq.10) then
      true(1)=0.2083229471647004E-05
      true(2)=0.8332935037760723E-11
      true(3)=0.9999979167621954E+00
      else if(i.eq.11) then
      true(1)=0.2083328471883087E-06
      true(2)=0.8333315602809495E-12
      true(3)=0.9999997916663195E+00
      else if(i.eq.12) then
      true(1)=0.2083340149701284E-07
      true(2)=0.8333360770334744E-13
      true(3)=0.9999999791665152E+00
      end if
      sum=0.0d+0
      do 270 k=1,3
      error(k)=dabs(true(k)-y(k))
      sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2)
270   continue
      sum=dsqrt(sum/dble(nd))
      sumh=sumh+sum
      summ=max(summ,sum)
 20     CONTINUE
CCC        CALL DTIME(TARRAY)
        it2=mclock()
        tarray(1)=(it2-it1)/100.0D+0
        WRITE(6,*)TARRAY(1),summ,sumh/12.0D+0
      summ=summ*atol
      sumh=sumh*atol/12.0D+0
      write(6,*) SUMM,SUMH
      us=log10(summ)
      ut=log10(sumh)
      write(6,*) us,ut
        ID=0
      NFE=IWORK(7)
      NJE=IWORK(8)
      NSTEP = IWORK(5)
      NDEC=IWORK(9)
      NBSOL=IWORK(10)
      NFAIL=IWORK(6)
        WRITE(6,901)NFE,NJE,NSTEP,NDEC,NBSOL,NFAIL
901   format(1x,6I12)
        WRITE(6,*)' ***** TOL=',RTOL,'  ELAPSED TIME=',TARRAY(1),' ****'
C -------- NEW TOLERANCE ---
 25     TOLST=TOLST*TOLFC
 30     CONTINUE
        STOP
        END

        SUBROUTINE F(N,X,Y,dy,IPAR,RPAR,ierr)
C --- RIGHT-HAND SIDE OF ROBERTSON EQUATION
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION Y(N),dy(N),IPAR(*),RPAR(*)
        dy(1)=-0.04D0*Y(1)+1.D4*Y(2)*Y(3)
        dy(3)=3.D7*Y(2)*Y(2)
        dy(2)=-dy(1)-dy(3)
        RETURN
        END 
C
        SUBROUTINE PDERV(X,Y,DFY,N,MEBND,IPAR,RPAR,ierr)
C --- JACOBIAN OF ROBERTSON EQUATION
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION Y(N),DFY(MEBND,N),IPAR(*),RPAR(*)
        PROD1=1.0D4*Y(2)
        PROD2=1.0D4*Y(3)
        PROD3=6.0D7*Y(2)
        DFY(1,1)=-0.04D0
        DFY(1,2)=PROD2
        DFY(1,3)=PROD1
        DFY(2,1)=0.04D0
        DFY(2,2)=-PROD2-PROD3
        DFY(2,3)=-PROD1
        DFY(3,1)=0.D0
        DFY(3,2)=PROD3
        DFY(3,3)=0.D0
        RETURN
        END

      Subroutine mas(n,am,m,IPAR,RPAR,ierr)
      dimension am(2,2),IPAR(*),RPAR(*)
      return
      end