C -----------------------------------------------------------------------------------
C     THIS MODULE IS PART OF THE CODE BIM.
C     THE CODE BIM NUMERICALLY SOLVES (STIFF) DIFFERENTIAL ODE 
C     PROBLEMS.
C
C     Copyright (C)2002-2007   
C
C     Authors: CECILIA MAGHERINI (cecilia.magherini@ing.unipi.it)
C              LUIGI   BRUGNANO  (brugnano@math.unifi.it) 
C
C
C     This program is free software; you can redistribute it and/or
C     modify it under the terms of the GNU General Public License
C     as published by the Free Software Foundation; either version 2
C     of the License, or (at your option) any later version.
C
C     This program is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     Licensed under The GNU General Public License, Version 2 or later.
C       http://www.gnu.org/licenses/info/GPLv2orLater.html
C
C     You should have received a copy of the GNU General Public License
C     along with this program; if not, write to the Free Software
C     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301,
C     USA.
C -----------------------------------------------------------------------------------


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                         c
c     In case the ISNAN function is not supported.        c
c                                                         c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      logical function isnan(A)
c      implicit none
c      double precision A, B
c      logical X
c      B = A 
c      X = ( A .EQ. B )
c      isnan = (.NOT.X)
c      return
c      end



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 
C     BLENDSTEP
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine blendstep4(m,y0,f0,Y,F,h,theta,ipvt,Z,
     &                      ldlu,mljac,mujac,ijob)
c
c     Blended iteration for the 4th order method
c      

      implicit none
      include 'param.dat'

      integer k
      parameter (k=3)


c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision  y0(m),f0(m),theta(ldlu,m),h

c     I/O parameters     
      double precision Y(m,k),F(m,k) 

c     Output parameters 
      double precision Z(m,k)

c     Local variables
      integer i,j

C ---------------------------------------------------------------------------------------
C 4th order BIM
C ---------------------------------------------------------------------------------------
      double precision DA4_1_1,DA4_1_2,DA4_1_3,DA4_1_4,
     &     DA4_2_1,DA4_2_2,DA4_2_3,DA4_2_4,
     &     DA4_3_1,DA4_3_2,DA4_3_3,DA4_3_4,
     &     DB4_1_1,DB4_1_2,DB4_1_3,DB4_1_4,
     &     DB4_2_1,DB4_2_2,DB4_2_3,DB4_2_4,
     &     DB4_3_1,DB4_3_2,DB4_3_3,DB4_3_4,
     &     A24_1_1,A24_1_2,A24_1_3,A24_1_4,
     &     A24_2_1,A24_2_2,A24_2_3,A24_2_4,
     &     A24_3_1,A24_3_2,A24_3_3,A24_3_4,
     &     B24_1_1,B24_2_1,B24_3_1

      parameter(
     & DA4_1_1=-102133D0/405000D0,  
     & DA4_1_2=98743D0/180000D0,  
     & DA4_1_3=-7387D0/22500D0,  
     & DA4_1_4=+51709D0/1620000D0,
     & DA4_2_1=-1D0-140353D0/810000D0,  
     & DA4_2_2=+7387D0/9000D0,     
     & DA4_2_3=10613D0/18000D0,   
     & DA4_2_4=-96031D0/405000D0, 
     & DA4_3_1=-22613D0/30000D0,   
     & DA4_3_2=-22161D0/20000D0, 
     & DA4_3_3=+22161D0/10000D0,   
     & DA4_3_4=-21257D0/60000D0,
     & A24_1_1=-302867D0/405000D0,  
     & A24_1_2= 81257D0/180000D0,  
     & A24_1_3= 7387D0/22500D0,  
     & A24_1_4=-51709D0/1620000D0,
     & A24_2_1= 140353D0/810000D0,  
     & A24_2_2=-7387D0/9000D0,    
     & A24_2_3= 7387D0/18000D0,   
     & A24_2_4= 96031D0/405000D0, 
     & A24_3_1=-7387D0/30000D0,   
     & A24_3_2= 22161D0/20000D0, 
     & A24_3_3=-22161D0/10000D0,   
     & A24_3_4= 81257D0/60000D0,
     & DB4_1_1= 919D0/13500D0,
     & DB4_1_2= 4589D0/30000D0,
     & DB4_1_3=-37D0/120D0,
     & DB4_1_4= 3D0/40D0, 
     & DB4_2_1= 115387D0/270000D0,
     & DB4_2_2= 17D0/15D0,
     & DB4_2_3=-6161D0/30000D0,
     & DB4_2_4=-1D0/15D0,
     & DB4_3_1= 3D0/8D0,
     & DB4_3_2= 9D0/8D0,
     & DB4_3_3= 9D0/8D0,
     & DB4_3_4=-3637D0/10000D0,
     & B24_1_1=7387D0/27000D0,
     & B24_2_1=-7387D0/270000D0,
     & B24_3_1=0D0)
C ---------------------------------------------------------------------------------------


c     Z=[y0 Y]*(DA)'
      do i=1,m
          Z(i,1)=y0(i)* DA4_1_1+Y(i,1)*DA4_1_2+Y(i,2)*DA4_1_3+
     &             Y(i,3)*DA4_1_4
            Z(i,2)=y0(i)* DA4_2_1+Y(i,1)*DA4_2_2+Y(i,2)*DA4_2_3+ 
     &             Y(i,3)*DA4_2_4
            Z(i,3)=y0(i)* DA4_3_1+Y(i,1)*DA4_3_2+Y(i,2)*DA4_3_3+
     &             Y(i,3)*DA4_3_4
      end do

c     Z=Z-h*[f0 F]*(dB)'
      do i=1,m 
           Z(i,1) = Z(i,1)-
     &           h*(f0(i)* DB4_1_1+F(i,1)*DB4_1_2+F(i,2)*DB4_1_3+
     &              F(i,3)*DB4_1_4)
           Z(i,2) = Z(i,2)-
     &           h*(f0(i)* DB4_2_1+F(i,1)*DB4_2_2+F(i,2)*DB4_2_3+
     &              F(i,3)*DB4_2_4)
           Z(i,3) = Z(i,3)-
     &           h*(f0(i)* DB4_3_1+F(i,1)*DB4_3_2+F(i,2)*DB4_3_3+
     &              F(i,3)*DB4_3_4)
      end do

c     theta*Z=Z
      call sollu(m,theta,ldlu,Z(1,1),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)

c     Z=Z+[y0 Y]*(A2)'
      do i=1,m
           Z(i,1)=Z(i,1)+ 
     &             y0(i)* A24_1_1+Y(i,1)*A24_1_2+Y(i,2)*A24_1_3+
     &              Y(i,3)* A24_1_4
             Z(i,2)=Z(i,2)+
     &               y0(i)* A24_2_1+Y(i,1)*A24_2_2+Y(i,2)*A24_2_3+ 
     &              Y(i,3)* A24_2_4
             Z(i,3)=Z(i,3)+
     &             y0(i)* A24_3_1+Y(i,1)*A24_3_2+Y(i,2)*A24_3_3+
     &              Y(i,3)* A24_3_4
      end do

c     Z=Z-h*[f0 F]*(B2)' 
      do i=1,m 
            Z(i,1)=Z(i,1)-h*(f0(i)*B24_1_1 + F(i,1)*gamma4)
            Z(i,2)=Z(i,2)-h*(f0(i)*B24_2_1 + F(i,2)*gamma4)
            Z(i,3)=Z(i,3)-h*(f0(i)*B24_3_1 + F(i,3)*gamma4)
      end do

c     theta*Z=Z
      call sollu(m,theta,ldlu,Z(1,1),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)

c
c     Y = Y-Z 
c
      do j=1,k
         do i=1,m
           Y(i,j) = Y(i,j) - Z(i,j)
         end do
      end do

      return
      end

C -------------------------------------------------------------------

      subroutine blendstep6(m,y0,f0,Y,F,h,theta,ipvt,Z,
     &                      ldlu,mljac,mujac,ijob)
c
c      Blended iteration for the 6th order method
c
      implicit none
      include 'param.dat'

      integer k
      parameter (k=4)


c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision y0(m),f0(m),theta(ldlu,m),h

c     I/O parameters
      double precision Y(m,k),F(m,k)

c     Output parameters
      double precision Z(m,k)

c     Local variables
      integer i,j

C ---------------------------------------------------------------------------------------
C 6th order BIM 
C ---------------------------------------------------------------------------------------

      double precision DA6_1_1,DA6_1_2,DA6_1_3,DA6_1_4,DA6_1_5,
     &     DA6_2_1,DA6_2_2,DA6_2_3,DA6_2_4,DA6_2_5,
     &     DA6_3_1,DA6_3_2,DA6_3_3,DA6_3_4,DA6_3_5,
     &     DA6_4_1,DA6_4_2,DA6_4_3,DA6_4_4,DA6_4_5,
     &     DB6_1_1,DB6_1_2,DB6_1_3,DB6_1_4,DB6_1_5,
     &     DB6_2_1,DB6_2_2,DB6_2_3,DB6_2_4,DB6_2_5,
     &     DB6_3_1,DB6_3_2,DB6_3_3,DB6_3_4,DB6_3_5,
     &     DB6_4_1,DB6_4_2,DB6_4_3,DB6_4_4,DB6_4_5,
     &     A26_1_1,A26_1_2,A26_1_3,A26_1_4,A26_1_5,
     &     A26_2_1,A26_2_2,A26_2_3,A26_2_4,A26_2_5,
     &     A26_3_1,A26_3_2,A26_3_3,A26_3_4,A26_3_5,
     &     A26_4_1,A26_4_2,A26_4_3,A26_4_4,A26_4_5,
     &     B26_1_1,B26_2_1,B26_3_1,B26_4_1

      parameter(
     & DA6_1_1=-1171629D0/5120000D0,
     & DA6_1_2=607997D0/960000D0,
     & DA6_1_3=-597981D0/1280000D0,
     & DA6_1_4=+4241D0/64000D0, 
     & DA6_1_5=-55133D0/15360000D0,
     & DA6_2_1=-307277D0/320000D0,
     & DA6_2_2=+4241D0/12000D0,
     & DA6_2_3=1D0+12723D0/80000D0,
     & DA6_2_4=-12723D0/20000D0,
     & DA6_2_5=+80579D0/960000D0,
     & DA6_3_1=-1D0-733693D0/5120000D0,
     & DA6_3_2=-4241D0/320000 D0,
     & DA6_3_3=+1234131D0/1280000D0,
     & DA6_3_4=137637D0/320000D0,
     & DA6_3_5=-1217167D0/5120000D0,
     & DA6_4_1=-24241D0/20000D0,
     & DA6_4_2=+4241D0/3750D0,
     & DA6_4_3=-12723D0/5000D0,
     & DA6_4_4=+4241D0/1250D0,
     & DA6_4_5=-1841D0/2400D0,
     & A26_1_1=-3948371D0/5120000D0,
     & A26_1_2=352003D0/960000D0,
     & A26_1_3=597981D0/1280000D0,
     & A26_1_4=-4241D0/64000D0, 
     & A26_1_5=55133D0/15360000D0,
     & A26_2_1=-12723D0/320000D0,
     & A26_2_2=-4241D0/12000D0,
     & A26_2_3=-12723D0/80000D0,
     & A26_2_4=12723D0/20000D0,
     & A26_2_5=-80579D0/960000D0,
     & A26_3_1=733693D0/5120000D0,
     & A26_3_2=4241D0/320000 D0,
     & A26_3_3=-1234131D0/1280000D0,
     & A26_3_4= 182363D0/320000D0,
     & A26_3_5= 1217167D0/5120000D0,
     & A26_4_1= 4241D0/20000D0,
     & A26_4_2=-4241D0/3750D0,
     & A26_4_3= 12723D0/5000D0,
     & A26_4_4=-4241D0/1250D0,
     & A26_4_5= 4241D0/2400D0,
     & DB6_1_1=60311D0/11520000D0,
     & DB6_1_2=7853D0/22500D0,
     & DB6_1_3=-49D0/60D0,
     & DB6_1_4=161D0/360D0,
     & DB6_1_5=-73D0/720D0,
     & DB6_2_1=257831D0/720000D0,
     & DB6_2_2=46D0/45D0,
     & DB6_2_3=-241D0/5000D0,
     & DB6_2_4=-14D0/45D0,
     & DB6_2_5=7D0/90D0,
     & DB6_3_1=1850413D0/3840000D0,
     & DB6_3_2=133D0/120D0,
     & DB6_3_3=23D0/20D0,
     & DB6_3_4=-1837D0/3750D0,
     & DB6_3_5=1D0/240D0,
     & DB6_4_1=14D0/45D0,
     & DB6_4_2=64D0/45D0,
     & DB6_4_3=8D0/15D0,
     & DB6_4_4=64D0/45D0, 
     & DB6_4_5=-24169D0/45000D0,
     & B26_1_1=343521D0/1280000D0,
     & B26_2_1=4241D0/80000D0,
     & B26_3_1=-131471D0/1280000D0,
     & B26_4_1=0D0)
C ---------------------------------------------------------------------------------------

c     Z=[y0 Y]*(DA)'
      do i=1,m
          Z(i,1)=y0(i)* DA6_1_1+Y(i,1)*DA6_1_2+Y(i,2)*DA6_1_3+
     &             Y(i,3)*DA6_1_4+Y(i,4)*DA6_1_5
            Z(i,2)=y0(i)* DA6_2_1+Y(i,1)*DA6_2_2+Y(i,2)*DA6_2_3+ 
     &             Y(i,3)*DA6_2_4+Y(i,4)*DA6_2_5
            Z(i,3)=y0(i)* DA6_3_1+Y(i,1)*DA6_3_2+Y(i,2)*DA6_3_3+
     &             Y(i,3)*DA6_3_4+Y(i,4)*DA6_3_5
            Z(i,4)=y0(i)* DA6_4_1+Y(i,1)*DA6_4_2+Y(i,2)*DA6_4_3+
     &             Y(i,3)*DA6_4_4+Y(i,4)*DA6_4_5
      end do

c     Z=Z-h*[f0 F]*(dB)'
      do i=1,m 
           Z(i,1) = Z(i,1)-
     &           h*(f0(i)* DB6_1_1+F(i,1)*DB6_1_2+F(i,2)*DB6_1_3+
     &              F(i,3)*DB6_1_4+F(i,4)*DB6_1_5)
           Z(i,2) = Z(i,2)-
     &           h*(f0(i)* DB6_2_1+F(i,1)*DB6_2_2+F(i,2)*DB6_2_3+
     &              F(i,3)*DB6_2_4+F(i,4)*DB6_2_5)
           Z(i,3) = Z(i,3)-
     &           h*(f0(i)* DB6_3_1+F(i,1)*DB6_3_2+F(i,2)*DB6_3_3+
     &              F(i,3)*DB6_3_4+F(i,4)*DB6_3_5)
           Z(i,4) = Z(i,4)-
     &           h*(f0(i)* DB6_4_1+F(i,1)*DB6_4_2+F(i,2)*DB6_4_3+
     &              F(i,3)*DB6_4_4+F(i,4)*DB6_4_5)
      end do

c     theta*Z=Z
      call sollu(m,theta,ldlu,Z(1,1),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,4),mljac,mujac,ipvt,ijob)

c     Z=Z+[y0 Y]*(A2)'
      do i=1,m
          Z(i,1)=Z(i,1)+
     &               y0(i)* A26_1_1+Y(i,1)*A26_1_2+Y(i,2)*A26_1_3+
     &              Y(i,3)* A26_1_4+Y(i,4)*A26_1_5
            Z(i,2)=Z(i,2)+
     &               y0(i)* A26_2_1+Y(i,1)*A26_2_2+Y(i,2)*A26_2_3+ 
     &              Y(i,3)* A26_2_4+Y(i,4)*A26_2_5
            Z(i,3)=Z(i,3)+
     &               y0(i)* A26_3_1+Y(i,1)*A26_3_2+Y(i,2)*A26_3_3+
     &              Y(i,3)* A26_3_4+Y(i,4)*A26_3_5
            Z(i,4)=Z(i,4)+
     &               y0(i)* A26_4_1+Y(i,1)*A26_4_2+Y(i,2)*A26_4_3+
     &              Y(i,3)* A26_4_4+Y(i,4)*A26_4_5
      end do

c     Z=Z-h*[f0 F]*(B2)' 
      do i=1,m 
            Z(i,1)=Z(i,1)-h*(f0(i)*B26_1_1 + F(i,1)*gamma6)
            Z(i,2)=Z(i,2)-h*(f0(i)*B26_2_1 + F(i,2)*gamma6)
            Z(i,3)=Z(i,3)-h*(f0(i)*B26_3_1 + F(i,3)*gamma6)
            Z(i,4)=Z(i,4)-h*(f0(i)*B26_4_1 + F(i,4)*gamma6)
      end do

c     theta*Z=Z
      call sollu(m,theta,ldlu,Z(1,1),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,4),mljac,mujac,ipvt,ijob)

c
c     Y = Y-Z 
c
      do j=1,k
         do i=1,m
           Y(i,j) = Y(i,j) - Z(i,j)
         end do
      end do

      return
      end

C -------------------------------------------------------------------

      subroutine blendstep8(m,y0,f0,Y,F,h,theta,ipvt,Z,
     &                      ldlu,mljac,mujac,ijob)
c
c      Blended iteration for the 8th order method
c
      implicit none
      include 'param.dat'
      integer k
      parameter (k=6)

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision y0(m),f0(m),theta(ldlu,m),h

c     I/O parameters
      double precision Y(m,k),F(m,k)

c     Output parameters
      double precision Z(m,k)

c     Local variables
      integer i,j

C ---------------------------------------------------------------------------------------
C 8th order BIM
C ---------------------------------------------------------------------------------------

      double precision 
     &     DA8_1_1,DA8_1_2,DA8_1_3,DA8_1_4,DA8_1_5,DA8_1_6,DA8_1_7,
     &     DA8_2_1,DA8_2_2,DA8_2_3,DA8_2_4,DA8_2_5,DA8_2_6,DA8_2_7,
     &     DA8_3_1,DA8_3_2,DA8_3_3,DA8_3_4,DA8_3_5,DA8_3_6,DA8_3_7,
     &     DA8_4_1,DA8_4_2,DA8_4_3,DA8_4_4,DA8_4_5,DA8_4_6,DA8_4_7,
     &     DA8_5_1,DA8_5_2,DA8_5_3,DA8_5_4,DA8_5_5,DA8_5_6,DA8_5_7,
     &     DA8_6_1,DA8_6_2,DA8_6_3,DA8_6_4,DA8_6_5,DA8_6_6,DA8_6_7,
     &     DB8_1_1,DB8_1_2,DB8_1_3,DB8_1_4,DB8_1_5,DB8_1_6,DB8_1_7,
     &     DB8_2_1,DB8_2_2,DB8_2_3,DB8_2_4,DB8_2_5,DB8_2_6,DB8_2_7,
     &     DB8_3_1,DB8_3_2,DB8_3_3,DB8_3_4,DB8_3_5,DB8_3_6,DB8_3_7,
     &     DB8_4_1,DB8_4_2,DB8_4_3,DB8_4_4,DB8_4_5,DB8_4_6,DB8_4_7,
     &     DB8_5_1,DB8_5_2,DB8_5_3,DB8_5_4,DB8_5_5,DB8_5_6,DB8_5_7,
     &     DB8_6_1,DB8_6_2,DB8_6_3,DB8_6_4,DB8_6_5,DB8_6_6,DB8_6_7,
     &     A28_1_1,A28_1_2,A28_1_3,A28_1_4,A28_1_5,A28_1_6,A28_1_7,
     &     A28_2_1,A28_2_2,A28_2_3,A28_2_4,A28_2_5,A28_2_6,A28_2_7,
     &     A28_3_1,A28_3_2,A28_3_3,A28_3_4,A28_3_5,A28_3_6,A28_3_7,
     &     A28_4_1,A28_4_2,A28_4_3,A28_4_4,A28_4_5,A28_4_6,A28_4_7,
     &     A28_5_1,A28_5_2,A28_5_3,A28_5_4,A28_5_5,A28_5_6,A28_5_7,
     &     A28_6_1,A28_6_2,A28_6_3,A28_6_4,A28_6_5,A28_6_6,A28_6_7,
     &     B28_1_1,B28_2_1,B28_3_1,B28_4_1,B28_5_1,B28_6_1

      parameter(
     &  DA8_1_1 =   -24312887D0/62208000D0,
     &  DA8_1_2 =   9595787D0/12960000D0,
     &  DA8_1_3 =      - 680419D0/2073600D0,
     &  DA8_1_4 =      - 263717D0/2332800D0,
     &  DA8_1_5 =      + 578429D0/4147200D0,
     &  DA8_1_6 =      - 147157D0/2592000D0,
     &  DA8_1_7 =      + 4150993D0/466560000D0,
     &  DA8_2_1 =   -1D0 - 394847D0/972000D0,
     &  DA8_2_2 =      +496837D0/405000D0,
     &  DA8_2_3 =   165733D0/648000D0,
     &  DA8_2_4 =      +24769D0/364500D0,
     &  DA8_2_5 =      -71393D0/324000D0,
     &  DA8_2_6 =      +1457D0/16200D0,
     &  DA8_2_7 =      -403589D0/29160000D0,
     &  DA8_3_1 =   -1D0 -33511D0/768000D0,
     &  DA8_3_2 =      +4371D0/160000D0,
     &  DA8_3_3 =      +48081D0/128000D0,
     &  DA8_3_4 =   1D0 +1457D0/9600D0,
     &  DA8_3_5 =      -161727D0/256000D0,
     &  DA8_3_6 =      +4371D0/32000D0,
     &  DA8_3_7 =      -10199D0/640000D0,
     &  DA8_4_1 =   -388381D0/486000D0,
     &  DA8_4_2 =      -85963D0/202500D0,
     &  DA8_4_3 =      +2914D0/10125D0,
     &  DA8_4_4 =      +71393D0/182250D0,
     &  DA8_4_5 =   145973D0/162000D0,
     &  DA8_4_6 =      -16027D0/40500D0,
     &  DA8_4_7 =      +141329D0/3645000D0,
     &  DA8_5_1 =   -1D0 -2073311D0/62208000D0,
     &  DA8_5_2 =      -106361D0/2592000D0,
     &  DA8_5_3 =      +893141D0/2073600D0,
     &  DA8_5_4 =      -2466701D0/2332800D0,
     &  DA8_5_5 =      +7187381D0/4147200D0,
     &  DA8_5_6 =   241859D0/2592000D0,
     &  DA8_5_7 =      -11695339D0/93312000D0,
     &  DA8_6_1 =   -1D0 -1457D0/12000D0,
     &  DA8_6_2 =      +4371D0/5000D0,
     &  DA8_6_3 =      -4371D0/1600D0,
     &  DA8_6_4 =      +1457D0/300D0,
     &  DA8_6_5 =      -4371D0/800D0,
     &  DA8_6_6 =      +4371D0/1000D0,
     &  DA8_6_7 =   -31393D0/40000D0,
     &  A28_1_1 =  -37895113D0/62208000D0,
     &  A28_1_2 =   3364213D0/12960000D0,
     &  A28_1_3 =   680419D0/2073600D0,
     &  A28_1_4 =   263717D0/2332800D0,
     &  A28_1_5 =   -578429D0/4147200D0,
     &  A28_1_6 =   147157D0/2592000D0,
     &  A28_1_7 =   -4150993D0/466560000D0,
     &  A28_2_1 =   394847D0/972000D0,
     &  A28_2_2 =   -496837D0/405000D0,
     &  A28_2_3 =   482267D0/648000D0,
     &  A28_2_4 =   -24769D0/364500D0,
     &  A28_2_5 =   71393D0/324000D0,
     &  A28_2_6 =   -1457D0/16200D0,
     &  A28_2_7 =   403589D0/29160000D0,
     &  A28_3_1 =   33511D0/768000D0,
     &  A28_3_2 =   -4371D0/160000D0,
     &  A28_3_3 =   -48081D0/128000D0,
     &  A28_3_4 =   -1457D0/9600D0,
     &  A28_3_5 =   161727D0/256000D0,
     &  A28_3_6 =   -4371D0/32000D0,
     &  A28_3_7 =   10199D0/640000D0,
     &  A28_4_1 =   -97619D0/486000D0,
     &  A28_4_2 =   85963D0/202500D0,
     &  A28_4_3 =   -2914D0/10125D0,
     &  A28_4_4 =   -71393D0/182250D0,
     &  A28_4_5 =   16027D0/162000D0,
     &  A28_4_6 =   16027D0/40500D0,
     &  A28_4_7 =   -141329D0/3645000D0,
     &  A28_5_1 =   2073311D0/62208000D0,
     &  A28_5_2 =   106361D0/2592000D0,
     &  A28_5_3 =   -893141D0/2073600D0,
     &  A28_5_4 =   2466701D0/2332800D0,
     &  A28_5_5 =   -7187381D0/4147200D0,
     &  A28_5_6 =   2350141D0/2592000D0,
     &  A28_5_7 =   11695339D0/93312000D0,
     &  A28_6_1 =   1457D0/12000D0,
     &  A28_6_2 =   -4371D0/5000D0,
     &  A28_6_3 =   4371D0/1600D0,
     &  A28_6_4 =   -1457D0/300D0,
     &  A28_6_5 =   4371D0/800D0,
     &  A28_6_6 =   -4371D0/1000D0,
     &  A28_6_7 =   71393D0/40000D0,
     &  DB8_1_1 =   87907D0/622080D0,
     &  DB8_1_2 =   3587D0/18000D0,
     &  DB8_1_3 =   -1141D0/2880D0,
     &  DB8_1_4 =   67D0/540D0,
     &  DB8_1_5 =   109D0/2880D0,
     &  DB8_1_6 =   -2D0/45D0,
     &  DB8_1_7 =   91D0/8640D0,
     &  DB8_2_1 =   203071D0/425250D0,
     &  DB8_2_2 =   2158D0/1575D0,
     &  DB8_2_3 =   -52291D0/126000D0,
     &  DB8_2_4 =   -52D0/945D0,
     &  DB8_2_5 =   23D0/252D0,
     &  DB8_2_6 =   -82D0/1575D0,
     &  DB8_2_7 =   199D0/18900D0,
     &  DB8_3_1 =   19177D0/64000D0,
     &  DB8_3_2 =   81D0/50D0,
     &  DB8_3_3 =   27D0/320D0,
     &  DB8_3_4 =   1643D0/2000D0,
     &  DB8_3_5 =   -243D0/320D0,
     &  DB8_3_6 =   27D0/100D0,
     &  DB8_3_7 =   -67D0/1600D0,
     &  DB8_4_1 =   716549D0/3402000D0,
     &  DB8_4_2 =   2368D0/1575D0,
     &  DB8_4_3 =   104D0/315D0,
     &  DB8_4_4 =   320D0/189D0,
     &  DB8_4_5 =   -78191D0/126000D0,
     &  DB8_4_6 =   128D0/1575D0,
     &  DB8_4_7 =   -64D0/4725D0,
     &  DB8_5_1 =   37053349D0/108864000D0,
     &  DB8_5_2 =   1739D0/1260D0,
     &  DB8_5_3 =   2713D0/4032D0,
     &  DB8_5_4 =   853D0/756D0,
     &  DB8_5_5 =   4463D0/4032D0,
     &  DB8_5_6 =   -40391D0/126000D0,
     &  DB8_5_7 =   -787D0/60480D0,
     &  DB8_6_1 =   41D0/140D0,
     &  DB8_6_2 =   54D0/35D0,
     &  DB8_6_3 =   27D0/140D0,
     &  DB8_6_4 =   68D0/35D0,
     &  DB8_6_5 =   27D0/140D0,
     &  DB8_6_6 =   54D0/35D0,
     &  DB8_6_7 =   -6099D0/14000D0,
     &  B28_1_1 =   24769D0/124416D0,
     &  B28_2_1 =   -18941D0/121500D0,
     &  B28_3_1 =   -1457D0/64000D0,
     &  B28_4_1 =   42253D0/486000D0,
     &  B28_5_1 =   -365707D0/15552000D0,
     &  B28_6_1 =   0D0)
C ---------------------------------------------------------------------------------------


c     Z=[y0 Y]*(DA)'
      do i=1,m
          Z(i,1)=y0(i)* DA8_1_1+Y(i,1)*DA8_1_2+Y(i,2)*DA8_1_3+
     &             Y(i,3)*DA8_1_4+Y(i,4)*DA8_1_5+Y(i,5)*DA8_1_6+
     &             Y(i,6)*DA8_1_7
            Z(i,2)=y0(i)* DA8_2_1+Y(i,1)*DA8_2_2+Y(i,2)*DA8_2_3+ 
     &             Y(i,3)*DA8_2_4+Y(i,4)*DA8_2_5+Y(i,5)*DA8_2_6+
     &             Y(i,6)*DA8_2_7
            Z(i,3)=y0(i)* DA8_3_1+Y(i,1)*DA8_3_2+Y(i,2)*DA8_3_3+
     &             Y(i,3)*DA8_3_4+Y(i,4)*DA8_3_5+Y(i,5)*DA8_3_6+
     &             Y(i,6)*DA8_3_7
            Z(i,4)=y0(i)* DA8_4_1+Y(i,1)*DA8_4_2+Y(i,2)*DA8_4_3+
     &             Y(i,3)*DA8_4_4+Y(i,4)*DA8_4_5+Y(i,5)*DA8_4_6+
     &             Y(i,6)*DA8_4_7
            Z(i,5)=y0(i)* DA8_5_1+Y(i,1)*DA8_5_2+Y(i,2)*DA8_5_3+
     &             Y(i,3)*DA8_5_4+Y(i,4)*DA8_5_5+Y(i,5)*DA8_5_6+
     &             Y(i,6)*DA8_5_7
            Z(i,6)=y0(i)* DA8_6_1+Y(i,1)*DA8_6_2+Y(i,2)*DA8_6_3+
     &             Y(i,3)*DA8_6_4+Y(i,4)*DA8_6_5+Y(i,5)*DA8_6_6+
     &             Y(i,6)*DA8_6_7
      end do

c     Z=Z-h*[f0 F]*(dB)'
      do i=1,m 
           Z(i,1) = Z(i,1)-
     &           h*(f0(i)* DB8_1_1+F(i,1)*DB8_1_2+F(i,2)*DB8_1_3+
     &              F(i,3)*DB8_1_4+F(i,4)*DB8_1_5+F(i,5)*DB8_1_6+
     &              F(i,6)*DB8_1_7)
           Z(i,2) = Z(i,2)-
     &           h*(f0(i)* DB8_2_1+F(i,1)*DB8_2_2+F(i,2)*DB8_2_3+
     &              F(i,3)*DB8_2_4+F(i,4)*DB8_2_5+F(i,5)*DB8_2_6+
     &              F(i,6)*DB8_2_7)
           Z(i,3) = Z(i,3)-
     &           h*(f0(i)* DB8_3_1+F(i,1)*DB8_3_2+F(i,2)*DB8_3_3+
     &              F(i,3)*DB8_3_4+F(i,4)*DB8_3_5+F(i,5)*DB8_3_6+
     &              F(i,6)*DB8_3_7)
           Z(i,4) = Z(i,4)-
     &           h*(f0(i)* DB8_4_1+F(i,1)*DB8_4_2+F(i,2)*DB8_4_3+
     &              F(i,3)*DB8_4_4+F(i,4)*DB8_4_5+F(i,5)*DB8_4_6+
     &              F(i,6)*DB8_4_7)
           Z(i,5) = Z(i,5)-
     &           h*(f0(i)* DB8_5_1+F(i,1)*DB8_5_2+F(i,2)*DB8_5_3+
     &              F(i,3)*DB8_5_4+F(i,4)*DB8_5_5+F(i,5)*DB8_5_6+
     &              F(i,6)*DB8_5_7)
           Z(i,6) = Z(i,6)-
     &           h*(f0(i)* DB8_6_1+F(i,1)*DB8_6_2+F(i,2)*DB8_6_3+
     &              F(i,3)*DB8_6_4+F(i,4)*DB8_6_5+F(i,5)*DB8_6_6+
     &              F(i,6)*DB8_6_7)
      end do

c     theta*Z=Z
      call sollu(m,theta,ldlu,Z(1,1),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,4),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,5),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,6),mljac,mujac,ipvt,ijob)

c     Z=Z+[y0 Y]*(A2)'
      do i=1,m
          Z(i,1)=Z(i,1)+
     &               y0(i)* A28_1_1+Y(i,1)*A28_1_2+Y(i,2)*A28_1_3+
     &              Y(i,3)* A28_1_4+Y(i,4)*A28_1_5+Y(i,5)*A28_1_6+
     &              Y(i,6)* A28_1_7 
            Z(i,2)=Z(i,2)+
     &               y0(i)* A28_2_1+Y(i,1)*A28_2_2+Y(i,2)*A28_2_3+ 
     &              Y(i,3)* A28_2_4+Y(i,4)*A28_2_5+Y(i,5)*A28_2_6+
     &              Y(i,6)* A28_2_7
            Z(i,3)=Z(i,3)+
     &               y0(i)* A28_3_1+Y(i,1)*A28_3_2+Y(i,2)*A28_3_3+
     &              Y(i,3)* A28_3_4+Y(i,4)*A28_3_5+Y(i,5)*A28_3_6+
     &              Y(i,6)* A28_3_7
            Z(i,4)=Z(i,4)+
     &               y0(i)* A28_4_1+Y(i,1)*A28_4_2+Y(i,2)*A28_4_3+
     &              Y(i,3)* A28_4_4+Y(i,4)*A28_4_5+Y(i,5)*A28_4_6+
     &              Y(i,6)* A28_4_7
            Z(i,5)=Z(i,5)+
     &               y0(i)* A28_5_1+Y(i,1)*A28_5_2+Y(i,2)*A28_5_3+
     &              Y(i,3)* A28_5_4+Y(i,4)*A28_5_5+Y(i,5)*A28_5_6+
     &              Y(i,6)* A28_5_7
            Z(i,6)=Z(i,6)+
     &               y0(i)* A28_6_1+Y(i,1)*A28_6_2+Y(i,2)*A28_6_3+
     &              Y(i,3)* A28_6_4+Y(i,4)*A28_6_5+Y(i,5)*A28_6_6+
     &              Y(i,6)* A28_6_7
      end do

c     Z=Z-h*[f0 F]*(B2)'
      do i=1,m 
            Z(i,1)=Z(i,1)-h*(f0(i)*B28_1_1 + F(i,1)*gamma8)
            Z(i,2)=Z(i,2)-h*(f0(i)*B28_2_1 + F(i,2)*gamma8)
            Z(i,3)=Z(i,3)-h*(f0(i)*B28_3_1 + F(i,3)*gamma8)
            Z(i,4)=Z(i,4)-h*(f0(i)*B28_4_1 + F(i,4)*gamma8)
            Z(i,5)=Z(i,5)-h*(f0(i)*B28_5_1 + F(i,5)*gamma8)
            Z(i,6)=Z(i,6)-h*(f0(i)*B28_6_1 + F(i,6)*gamma8)
      end do

c     theta*Z=Z
      call sollu(m,theta,ldlu,Z(1,1),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,4),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,5),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,6),mljac,mujac,ipvt,ijob)

c
c     Y = Y-Z 
c
      do j=1,k
         do i=1,m
           Y(i,j) = Y(i,j) - Z(i,j)
         end do
      end do


      return
      end

C -------------------------------------------------------------------

      subroutine blendstep10(m,y0,f0,Y,F,h,theta,ipvt,Z,
     &                      ldlu,mljac,mujac,ijob)
c
c      Blended iteration for the 10th order method
c
      implicit none
      include 'param.dat'
      integer k
      parameter (k=8)

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision y0(m),f0(m),theta(ldlu,m),h

c     I/O parameters
      double precision Y(m,k),F(m,k)

c     Output parameters
      double precision Z(m,k)

c     Local variables
      integer i,j

C ---------------------------------------------------------------------------------------
C 10th order BIM
C ---------------------------------------------------------------------------------------

      double precision 
     &   DA10_1_1,DA10_1_2,DA10_1_3,DA10_1_4,DA10_1_5,DA10_1_6,DA10_1_7,
     &   DA10_1_8,DA10_1_9,
     &   DA10_2_1,DA10_2_2,DA10_2_3,DA10_2_4,DA10_2_5,DA10_2_6,DA10_2_7,
     &   DA10_2_8,DA10_2_9,
     &   DA10_3_1,DA10_3_2,DA10_3_3,DA10_3_4,DA10_3_5,DA10_3_6,DA10_3_7,
     &   DA10_3_8,DA10_3_9,
     &   DA10_4_1,DA10_4_2,DA10_4_3,DA10_4_4,DA10_4_5,DA10_4_6,DA10_4_7,
     &   DA10_4_8,DA10_4_9,
     &   DA10_5_1,DA10_5_2,DA10_5_3,DA10_5_4,DA10_5_5,DA10_5_6,DA10_5_7,
     &   DA10_5_8,DA10_5_9,
     &   DA10_6_1,DA10_6_2,DA10_6_3,DA10_6_4,DA10_6_5,DA10_6_6,DA10_6_7,
     &   DA10_6_8,DA10_6_9,
     &   DA10_7_1,DA10_7_2,DA10_7_3,DA10_7_4,DA10_7_5,DA10_7_6,DA10_7_7,
     &   DA10_7_8,DA10_7_9,
     &   DA10_8_1,DA10_8_2,DA10_8_3,DA10_8_4,DA10_8_5,DA10_8_6,DA10_8_7,
     &   DA10_8_8,DA10_8_9,
     &   DB10_1_1,DB10_1_2,DB10_1_3,DB10_1_4,DB10_1_5,DB10_1_6,DB10_1_7,
     &   DB10_1_8,DB10_1_9,
     &   DB10_2_1,DB10_2_2,DB10_2_3,DB10_2_4,DB10_2_5,DB10_2_6,DB10_2_7,
     &   DB10_2_8,DB10_2_9,
     &   DB10_3_1,DB10_3_2,DB10_3_3,DB10_3_4,DB10_3_5,DB10_3_6,DB10_3_7,
     &   DB10_3_8,DB10_3_9,
     &   DB10_4_1,DB10_4_2,DB10_4_3,DB10_4_4,DB10_4_5,DB10_4_6,DB10_4_7,
     &   DB10_4_8,DB10_4_9,
     &   DB10_5_1,DB10_5_2,DB10_5_3,DB10_5_4,DB10_5_5,DB10_5_6,DB10_5_7,
     &   DB10_5_8,DB10_5_9,
     &   DB10_6_1,DB10_6_2,DB10_6_3,DB10_6_4,DB10_6_5,DB10_6_6,DB10_6_7,
     &   DB10_6_8,DB10_6_9,
     &   DB10_7_1,DB10_7_2,DB10_7_3,DB10_7_4,DB10_7_5,DB10_7_6,DB10_7_7,
     &   DB10_7_8,DB10_7_9,
     &   DB10_8_1,DB10_8_2,DB10_8_3,DB10_8_4,DB10_8_5,DB10_8_6,DB10_8_7,
     &   DB10_8_8,DB10_8_9,
     &   A210_1_1,A210_1_2,A210_1_3,A210_1_4,A210_1_5,A210_1_6,A210_1_7,
     &   A210_1_8,A210_1_9,
     &   A210_2_1,A210_2_2,A210_2_3,A210_2_4,A210_2_5,A210_2_6,A210_2_7,
     &   A210_2_8,A210_2_9,
     &   A210_3_1,A210_3_2,A210_3_3,A210_3_4,A210_3_5,A210_3_6,A210_3_7,
     &   A210_3_8,A210_3_9,
     &   A210_4_1,A210_4_2,A210_4_3,A210_4_4,A210_4_5,A210_4_6,A210_4_7,
     &   A210_4_8,A210_4_9,
     &   A210_5_1,A210_5_2,A210_5_3,A210_5_4,A210_5_5,A210_5_6,A210_5_7,
     &   A210_5_8,A210_5_9,
     &   A210_6_1,A210_6_2,A210_6_3,A210_6_4,A210_6_5,A210_6_6,A210_6_7,
     &   A210_6_8,A210_6_9,
     &   A210_7_1,A210_7_2,A210_7_3,A210_7_4,A210_7_5,A210_7_6,A210_7_7,
     &   A210_7_8,A210_7_9,
     &   A210_8_1,A210_8_2,A210_8_3,A210_8_4,A210_8_5,A210_8_6,A210_8_7,
     &   A210_8_8,A210_8_9,
     &   B210_1_1,B210_2_1,B210_3_1,B210_4_1,B210_5_1,B210_6_1,B210_7_1,
     &   B210_8_1

 
      parameter(
     &  DA10_1_1    =   -1116185640857D0/1342177280000D0,
     &  DA10_1_2    =   1D0 +121395173141D0/146800640000D0,
     &  DA10_1_3    =   -6468370013D0/3355443200D0,
     &  DA10_1_4    =   +22440146897D0/12582912000D0,
     &  DA10_1_5    =   -57418737523D0/40265318400D0,
     &  DA10_1_6    =   +17489295313D0/20971520000D0,
     &  DA10_1_7    =   -82495624981D0/251658240000D0,
     &  DA10_1_8    =   +193947079D0/2516582400D0,
     &  DA10_1_9    =   -15367501777D0/1879048192000D0,
     &  DA10_2_1    =   -1D0    -359763461D0/1048576000D0,
     &  DA10_2_2    =   +133859921D0/114688000D0,
     &  DA10_2_3    =   -21140373D0/327680000D0,
     &  DA10_2_4    =   +15150619D0/16384000D0,
     &  DA10_2_5    =   -67553873D0/52428800D0,
     &  DA10_2_6    =   +44952727D0/49152000D0,
     &  DA10_2_7    =   -26201627D0/65536000D0,
     &  DA10_2_8    =   +8246437D0/81920000D0,
     &  DA10_2_9    =   -246865651D0/22020096000D0,
     &  DA10_3_1    =   -552495440623D0/805306368000D0,
     &  DA10_3_2    =   -28197945999D0/29360128000D0,
     &  DA10_3_3    =   +32441533071D0/16777216000D0,
     &  DA10_3_4    =   -51906383051D0/62914560000D0,
     &  DA10_3_5    =   +15467184783D0/13421772800D0,
     &  DA10_3_6    =   -3942154371D0/4194304000D0,
     &  DA10_3_7    =   +21125300879D0/50331648000D0,
     &  DA10_3_8    =   -445562559D0/4194304000D0,
     &  DA10_3_9    =   +111284620491D0/9395240960000D0,
     &  DA10_4_1    =   -19796057D0/20480000D0,
     &  DA10_4_2    =   -535553D0/6720000D0,
     &  DA10_4_3    =   +63403D0/1280000D0,
     &  DA10_4_4    =   +281941D0/960000D0,
     &  DA10_4_5    =   1D0 +9443D0/40960D0,
     &  DA10_4_6    =   -219887D0/320000D0,
     &  DA10_4_7    =   +754091D0/3840000D0,
     &  DA10_4_8    =   -39121D0/960000D0,
     &  DA10_4_9    =   +581419D0/143360000D0,
     &  DA10_5_1    =   -1D0    -10781441377D0/53687091200D0,
     &  DA10_5_2    =   +16856063921D0/29360128000D0,
     &  DA10_5_3    =   -47031869203D0/50331648000D0,
     &  DA10_5_4    =   +4439913787D0/4194304000D0,
     &  DA10_5_5    =   -6247096241D0/13421772800D0,
     &  DA10_5_6    =   1D0 +447917913D0/838860800D0,
     &  DA10_5_7    =   -11512558907D0/16777216000D0,
     &  DA10_5_8    =   +560510849D0/4194304000D0,
     &  DA10_5_9    =   -75322449683D0/5637144576000D0,
     &  DA10_6_1    =   -14564721451D0/15728640000D0,
     &  DA10_6_2    =   -109580619D0/573440000D0,
     &  DA10_6_3    =   +15285519D0/65536000D0,
     &  DA10_6_4    =   -4234511D0/49152000D0,
     &  DA10_6_5    =   -3573501D0/10485760D0,
     &  DA10_6_6    =   +84173553D0/81920000D0,
     &  DA10_6_7    =   484815179D0/983040000D0,
     &  DA10_6_8    =   -3694911D0/16384000D0,
     &  DA10_6_9    =   +114744591D0/7340032000D0,
     &  DA10_7_1    =   -10378315789D0/10737418240D0,
     &  DA10_7_2    =   -2207113739D0/12582912000D0,
     &  DA10_7_3    =   +48853634347D0/83886080000D0,
     &  DA10_7_4    =   -16701919087D0/12582912000D0,
     &  DA10_7_5    =   +84761718349D0/40265318400D0,
     &  DA10_7_6    =   -2054315207D0/838860800D0,
     &  DA10_7_7    =   +120669802351D0/50331648000D0,
     &  DA10_7_8    =   -1748617243D0/20971520000D0,
     &  DA10_7_9    =   -22368269479D0/268435456000D0,
     &  DA10_8_1    =   -1D0    -1349D0/16000D0,
     &  DA10_8_2    =   +1349D0/1750D0,
     &  DA10_8_3    =   -9443D0/3000D0,
     &  DA10_8_4    =   +9443D0/1250D0,
     &  DA10_8_5    =   -9443D0/800D0,
     &  DA10_8_6    =   +9443D0/750D0,
     &  DA10_8_7    =   -9443D0/1000D0,
     &  DA10_8_8    =   +1349D0/250D0,
     &  DA10_8_9    =   -466589D0/560000D0)
        parameter(
     &  A210_1_1    =   -225991639143D0/1342177280000D0,
     &  A210_1_2    =   -121395173141D0/146800640000D0,
     &  A210_1_3    =   +6468370013D0/3355443200D0,
     &  A210_1_4    =   -22440146897D0/12582912000D0,
     &  A210_1_5    =   +57418737523D0/40265318400D0,
     &  A210_1_6    =   -17489295313D0/20971520000D0,
     &  A210_1_7    =   +82495624981D0/251658240000D0,
     &  A210_1_8    =   -193947079D0/2516582400D0,
     &  A210_1_9    =   +15367501777D0/1879048192000D0,
     &  A210_2_1    =   +359763461D0/1048576000D0,
     &  A210_2_2    =   -133859921D0/114688000D0,
     &  A210_2_3    =   +348820373D0/327680000D0,
     &  A210_2_4    =   -15150619D0/16384000D0,
     &  A210_2_5    =   +67553873D0/52428800D0,
     &  A210_2_6    =   -44952727D0/49152000D0,
     &  A210_2_7    =   +26201627D0/65536000D0,
     &  A210_2_8    =   -8246437D0/81920000D0,
     &  A210_2_9    =   +246865651D0/22020096000D0,
     &  A210_3_1    =   -252810927377D0/805306368000D0,
     &  A210_3_2    =   +28197945999D0/29360128000D0,
     &  A210_3_3    =   -32441533071D0/16777216000D0,
     &  A210_3_4    =   +114820943051D0/62914560000D0,
     &  A210_3_5    =   -15467184783D0/13421772800D0,
     &  A210_3_6    =   +3942154371D0/4194304000D0,
     &  A210_3_7    =   -21125300879D0/50331648000D0,
     &  A210_3_8    =   +445562559D0/4194304000D0,
     &  A210_3_9    =   -111284620491D0/9395240960000D0,
     &  A210_4_1    =   -683943D0/20480000D0,
     &  A210_4_2    =   +535553D0/6720000D0,
     &  A210_4_3    =   -63403D0/1280000D0,
     &  A210_4_4    =   -281941D0/960000D0,
     &  A210_4_5    =   -9443D0/40960D0,
     &  A210_4_6    =   +219887D0/320000D0,
     &  A210_4_7    =   -754091D0/3840000D0,
     &  A210_4_8    =   +39121D0/960000D0,
     &  A210_4_9    =   -581419D0/143360000D0,
     &  A210_5_1    =   +10781441377D0/53687091200D0,
     &  A210_5_2    =   -16856063921D0/29360128000D0,
     &  A210_5_3    =   +47031869203D0/50331648000D0,
     &  A210_5_4    =   -4439913787D0/4194304000D0,
     &  A210_5_5    =   +6247096241D0/13421772800D0,
     &  A210_5_6    =   -447917913D0/838860800D0,
     &  A210_5_7    =   +11512558907D0/16777216000D0,
     &  A210_5_8    =   -560510849D0/4194304000D0,
     &  A210_5_9    =   +75322449683D0/5637144576000D0,
     &  A210_6_1    =   -1163918549D0/15728640000D0,
     &  A210_6_2    =   +109580619D0/573440000D0,
     &  A210_6_3    =   -15285519D0/65536000D0,
     &  A210_6_4    =   +4234511D0/49152000D0,
     &  A210_6_5    =   +3573501D0/10485760D0,
     &  A210_6_6    =   -84173553D0/81920000D0,
     &  A210_6_7    =   +498224821D0/983040000D0,
     &  A210_6_8    =   +3694911D0/16384000D0,
     &  A210_6_9    =   -114744591D0/7340032000D0,
     &  A210_7_1    =   -359102451D0/10737418240D0,
     &  A210_7_2    =   +2207113739D0/12582912000D0,
     &  A210_7_3    =   -48853634347D0/83886080000D0,
     &  A210_7_4    =   +16701919087D0/12582912000D0,
     &  A210_7_5    =   -84761718349D0/40265318400D0,
     &  A210_7_6    =   +2054315207D0/838860800D0,
     &  A210_7_7    =   -120669802351D0/50331648000D0,
     &  A210_7_8    =   +22720137243D0/20971520000D0,
     &  A210_7_9    =   +22368269479D0/268435456000D0,
     &  A210_8_1    =   +1349D0/16000D0,
     &  A210_8_2    =   -1349D0/1750D0,
     &  A210_8_3    =   +9443D0/3000D0,
     &  A210_8_4    =   -9443D0/1250D0,
     &  A210_8_5    =   +9443D0/800D0,
     &  A210_8_6    =   -9443D0/750D0,
     &  A210_8_7    =   +9443D0/1000D0,
     &  A210_8_8    =   -1349D0/250D0)
        parameter(
     &  A210_8_9    =   +1026589D0/560000D0,
     &  DB10_1_1    =   +761073292818799D0/2720626900992000D0,
     &  DB10_1_2    =   +557937353D0/1297296000D0,
     &  DB10_1_3    =   -214412651D0/259459200D0,
     &  DB10_1_4    =   +170436457D0/259459200D0,
     &  DB10_1_5    =   -7282853D0/25945920D0,
     &  DB10_1_6    =   -4668473D0/259459200D0,
     &  DB10_1_7    =   +22494019D0/259459200D0,
     &  DB10_1_8    =   -10434029D0/259459200D0,
     &  DB10_1_9    =   +3345851D0/518918400D0,
     &  DB10_2_1    =   +1183625282033D0/2975685672960D0,
     &  DB10_2_2    =   +23930143D0/14189175D0,
     &  DB10_2_3    =   -213517489D0/162162000D0,
     &  DB10_2_4    =   +3230893D0/2027025D0,
     &  DB10_2_5    =   -686209D0/405405D0,
     &  DB10_2_6    =   +2420083D0/2027025D0,
     &  DB10_2_7    =   -4396151D0/8108100D0,
     &  DB10_2_8    =   +2038273D0/14189175D0,
     &  DB10_2_9    =   -1918153D0/113513400D0,
     &  DB10_3_1    =   +70530807725549D0/423208629043200D0,
     &  DB10_3_2    =   +330413059D0/201801600D0,
     &  DB10_3_3    =   -3305147D0/28828800D0,
     &  DB10_3_4    =   +30733931D0/20592000D0,
     &  DB10_3_5    =   -4993271D0/2882880D0,
     &  DB10_3_6    =   +4741897D0/4118400D0,
     &  DB10_3_7    =   -14566397D0/28828800D0,
     &  DB10_3_8    =   +26359309D0/201801600D0,
     &  DB10_3_9    =   -6083071D0/403603200D0,
     &  DB10_4_1    =   +3304437553D0/11623772160D0,
     &  DB10_4_2    =   +21369776D0/14189175D0,
     &  DB10_4_3    =   +92756D0/289575D0,
     &  DB10_4_4    =   +3454736D0/2027025D0,
     &  DB10_4_5    =   -90373469D0/162162000D0,
     &  DB10_4_6    =   +91376D0/2027025D0,
     &  DB10_4_7    =   +6956D0/289575D0,
     &  DB10_4_8    =   -251824D0/14189175D0,
     &  DB10_4_9    =   +46493D0/14189175D0,
     &  DB10_5_1    =   +6901687296484073D0/19044388306944000D0,
     &  DB10_5_2    =   +115012739D0/72648576D0,
     &  DB10_5_3    =   +628493D0/10378368D0,
     &  DB10_5_4    =   +22584689D0/10378368D0,
     &  DB10_5_5    =   -4595D0/741312D0,
     &  DB10_5_6    =   +630493723D0/1297296000D0,
     &  DB10_5_7    =   -3425557D0/10378368D0,
     &  DB10_5_8    =   +5553389D0/72648576D0,
     &  DB10_5_9    =   -1181891D0/145297152D0,
     &  DB10_6_1    =   +160618157669D0/635830272000D0,
     &  DB10_6_2    =   +198307D0/121275D0,
     &  DB10_6_3    =   -9029D0/69300D0,
     &  DB10_6_4    =   +44857D0/17325D0,
     &  DB10_6_5    =   -2131D0/3465D0,
     &  DB10_6_6    =   +37927D0/17325D0,
     &  DB10_6_7    =   -976837D0/1386000D0,
     &  DB10_6_8    =   +11197D0/121275D0,
     &  DB10_6_9    =   -11197D0/970200D0,
     &  DB10_7_1    =   +153681594740771D0/544125380198400D0,
     &  DB10_7_2    =   +404839597D0/259459200D0,
     &  DB10_7_3    =   +5425339D0/37065600D0,
     &  DB10_7_4    =   +73758847D0/37065600D0,
     &  DB10_7_5    =   +878107D0/3706560D0,
     &  DB10_7_6    =   +48743857D0/37065600D0,
     &  DB10_7_7    =   +39269149D0/37065600D0,
     &  DB10_7_8    =   -335961817D0/1297296000D0,
     &  DB10_7_9    =   -7219753D0/518918400D0,
     &  DB10_8_1    =   +3956D0/14175D0,
     &  DB10_8_2    =   +23552D0/14175D0,
     &  DB10_8_3    =   -3712D0/14175D0,
     &  DB10_8_4    =   +41984D0/14175D0,
     &  DB10_8_5    =   -3632D0/2835D0,
     &  DB10_8_6    =   +41984D0/14175D0,
     &  DB10_8_7    =   -3712D0/14175D0,
     &  DB10_8_8    =   +23552D0/14175D0,
     &  DB10_8_9    =   -448403D0/1134000D0,
     &  B210_1_1    =   +1037851801D0/33554432000D0,
     &  B210_2_1    =   -3193083D0/26214400D0,
     &  B210_3_1    =   +153048097D0/1342177280D0,
     &  B210_4_1    =   +1349D0/102400D0,
     &  B210_5_1    =   -2509046919D0/33554432000D0,
     &  B210_6_1    =   +3762361D0/131072000D0,
     &  B210_7_1    =   +52838981D0/6710886400D0,
     &  B210_8_1    =   0D0)
C ---------------------------------------------------------------------------------------

c     Z=[y0 Y]*(DA)'
      do i=1,m
          Z(i,1)=y0(i)* DA10_1_1+Y(i,1)*DA10_1_2+Y(i,2)*DA10_1_3+
     &             Y(i,3)*DA10_1_4+Y(i,4)*DA10_1_5+Y(i,5)*DA10_1_6+
     &             Y(i,6)*DA10_1_7+Y(i,7)*DA10_1_8+Y(i,8)*DA10_1_9
            Z(i,2)=y0(i)* DA10_2_1+Y(i,1)*DA10_2_2+Y(i,2)*DA10_2_3+ 
     &             Y(i,3)*DA10_2_4+Y(i,4)*DA10_2_5+Y(i,5)*DA10_2_6+
     &             Y(i,6)*DA10_2_7+Y(i,7)*DA10_2_8+Y(i,8)*DA10_2_9
            Z(i,3)=y0(i)* DA10_3_1+Y(i,1)*DA10_3_2+Y(i,2)*DA10_3_3+
     &             Y(i,3)*DA10_3_4+Y(i,4)*DA10_3_5+Y(i,5)*DA10_3_6+
     &             Y(i,6)*DA10_3_7+Y(i,7)*DA10_3_8+Y(i,8)*DA10_3_9
            Z(i,4)=y0(i)* DA10_4_1+Y(i,1)*DA10_4_2+Y(i,2)*DA10_4_3+
     &             Y(i,3)*DA10_4_4+Y(i,4)*DA10_4_5+Y(i,5)*DA10_4_6+
     &             Y(i,6)*DA10_4_7+Y(i,7)*DA10_4_8+Y(i,8)*DA10_4_9
            Z(i,5)=y0(i)* DA10_5_1+Y(i,1)*DA10_5_2+Y(i,2)*DA10_5_3+
     &             Y(i,3)*DA10_5_4+Y(i,4)*DA10_5_5+Y(i,5)*DA10_5_6+
     &             Y(i,6)*DA10_5_7+Y(i,7)*DA10_5_8+Y(i,8)*DA10_5_9
            Z(i,6)=y0(i)* DA10_6_1+Y(i,1)*DA10_6_2+Y(i,2)*DA10_6_3+
     &             Y(i,3)*DA10_6_4+Y(i,4)*DA10_6_5+Y(i,5)*DA10_6_6+
     &             Y(i,6)*DA10_6_7+Y(i,7)*DA10_6_8+Y(i,8)*DA10_6_9
            Z(i,7)=y0(i)* DA10_7_1+Y(i,1)*DA10_7_2+Y(i,2)*DA10_7_3+
     &             Y(i,3)*DA10_7_4+Y(i,4)*DA10_7_5+Y(i,5)*DA10_7_6+
     &             Y(i,6)*DA10_7_7+Y(i,7)*DA10_7_8+Y(i,8)*DA10_7_9
            Z(i,8)=y0(i)* DA10_8_1+Y(i,1)*DA10_8_2+Y(i,2)*DA10_8_3+
     &             Y(i,3)*DA10_8_4+Y(i,4)*DA10_8_5+Y(i,5)*DA10_8_6+
     &             Y(i,6)*DA10_8_7+Y(i,7)*DA10_8_8+Y(i,8)*DA10_8_9
      end do

c     Z=Z-h*[f0 F]*(dB)'
      do i=1,m 
           Z(i,1) = Z(i,1)-
     &           h*(f0(i)* DB10_1_1+F(i,1)*DB10_1_2+F(i,2)*DB10_1_3+
     &              F(i,3)*DB10_1_4+F(i,4)*DB10_1_5+F(i,5)*DB10_1_6+
     &              F(i,6)*DB10_1_7+F(i,7)*DB10_1_8+F(i,8)*DB10_1_9)
           Z(i,2) = Z(i,2)-
     &           h*(f0(i)* DB10_2_1+F(i,1)*DB10_2_2+F(i,2)*DB10_2_3+
     &              F(i,3)*DB10_2_4+F(i,4)*DB10_2_5+F(i,5)*DB10_2_6+
     &              F(i,6)*DB10_2_7+F(i,7)*DB10_2_8+F(i,8)*DB10_2_9)
           Z(i,3) = Z(i,3)-
     &           h*(f0(i)* DB10_3_1+F(i,1)*DB10_3_2+F(i,2)*DB10_3_3+
     &              F(i,3)*DB10_3_4+F(i,4)*DB10_3_5+F(i,5)*DB10_3_6+
     &              F(i,6)*DB10_3_7+F(i,7)*DB10_3_8+F(i,8)*DB10_3_9)
           Z(i,4) = Z(i,4)-
     &           h*(f0(i)* DB10_4_1+F(i,1)*DB10_4_2+F(i,2)*DB10_4_3+
     &              F(i,3)*DB10_4_4+F(i,4)*DB10_4_5+F(i,5)*DB10_4_6+
     &              F(i,6)*DB10_4_7+F(i,7)*DB10_4_8+F(i,8)*DB10_4_9)
           Z(i,5) = Z(i,5)-
     &           h*(f0(i)* DB10_5_1+F(i,1)*DB10_5_2+F(i,2)*DB10_5_3+
     &              F(i,3)*DB10_5_4+F(i,4)*DB10_5_5+F(i,5)*DB10_5_6+
     &              F(i,6)*DB10_5_7+F(i,7)*DB10_5_8+F(i,8)*DB10_5_9)
           Z(i,6) = Z(i,6)-
     &           h*(f0(i)* DB10_6_1+F(i,1)*DB10_6_2+F(i,2)*DB10_6_3+
     &              F(i,3)*DB10_6_4+F(i,4)*DB10_6_5+F(i,5)*DB10_6_6+
     &              F(i,6)*DB10_6_7+F(i,7)*DB10_6_8+F(i,8)*DB10_6_9)
           Z(i,7) = Z(i,7)-
     &           h*(f0(i)* DB10_7_1+F(i,1)*DB10_7_2+F(i,2)*DB10_7_3+
     &              F(i,3)*DB10_7_4+F(i,4)*DB10_7_5+F(i,5)*DB10_7_6+
     &              F(i,6)*DB10_7_7+F(i,7)*DB10_7_8+F(i,8)*DB10_7_9)
           Z(i,8) = Z(i,8)-
     &           h*(f0(i)* DB10_8_1+F(i,1)*DB10_8_2+F(i,2)*DB10_8_3+
     &              F(i,3)*DB10_8_4+F(i,4)*DB10_8_5+F(i,5)*DB10_8_6+
     &              F(i,6)*DB10_8_7+F(i,7)*DB10_8_8+F(i,8)*DB10_8_9)
      end do

c     theta*Z=Z
      call sollu(m,theta,ldlu,Z(1,1),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,4),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,5),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,6),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,7),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,8),mljac,mujac,ipvt,ijob)

c     Z=Z+[y0 Y]*(A2)'
      do i=1,m
          Z(i,1)=Z(i,1)+
     &               y0(i)* A210_1_1+Y(i,1)*A210_1_2+Y(i,2)*A210_1_3+
     &              Y(i,3)* A210_1_4+Y(i,4)*A210_1_5+Y(i,5)*A210_1_6+
     &              Y(i,6)* A210_1_7+Y(i,7)*A210_1_8+Y(i,8)*A210_1_9
            Z(i,2)=Z(i,2)+
     &               y0(i)* A210_2_1+Y(i,1)*A210_2_2+Y(i,2)*A210_2_3+ 
     &              Y(i,3)* A210_2_4+Y(i,4)*A210_2_5+Y(i,5)*A210_2_6+
     &              Y(i,6)* A210_2_7+Y(i,7)*A210_2_8+Y(i,8)*A210_2_9
            Z(i,3)=Z(i,3)+
     &               y0(i)* A210_3_1+Y(i,1)*A210_3_2+Y(i,2)*A210_3_3+
     &              Y(i,3)* A210_3_4+Y(i,4)*A210_3_5+Y(i,5)*A210_3_6+
     &              Y(i,6)* A210_3_7+Y(i,7)*A210_3_8+Y(i,8)*A210_3_9
            Z(i,4)=Z(i,4)+
     &               y0(i)* A210_4_1+Y(i,1)*A210_4_2+Y(i,2)*A210_4_3+
     &              Y(i,3)* A210_4_4+Y(i,4)*A210_4_5+Y(i,5)*A210_4_6+
     &              Y(i,6)* A210_4_7+Y(i,7)*A210_4_8+Y(i,8)*A210_4_9
            Z(i,5)=Z(i,5)+
     &               y0(i)* A210_5_1+Y(i,1)*A210_5_2+Y(i,2)*A210_5_3+
     &              Y(i,3)* A210_5_4+Y(i,4)*A210_5_5+Y(i,5)*A210_5_6+
     &              Y(i,6)* A210_5_7+Y(i,7)*A210_5_8+Y(i,8)*A210_5_9
            Z(i,6)=Z(i,6)+
     &               y0(i)* A210_6_1+Y(i,1)*A210_6_2+Y(i,2)*A210_6_3+
     &              Y(i,3)* A210_6_4+Y(i,4)*A210_6_5+Y(i,5)*A210_6_6+
     &              Y(i,6)* A210_6_7+Y(i,7)*A210_6_8+Y(i,8)*A210_6_9
            Z(i,7)=Z(i,7)+
     &               y0(i)* A210_7_1+Y(i,1)*A210_7_2+Y(i,2)*A210_7_3+
     &              Y(i,3)* A210_7_4+Y(i,4)*A210_7_5+Y(i,5)*A210_7_6+
     &              Y(i,6)* A210_7_7+Y(i,7)*A210_7_8+Y(i,8)*A210_7_9
            Z(i,8)=Z(i,8)+
     &               y0(i)* A210_8_1+Y(i,1)*A210_8_2+Y(i,2)*A210_8_3+
     &              Y(i,3)* A210_8_4+Y(i,4)*A210_8_5+Y(i,5)*A210_8_6+
     &              Y(i,6)* A210_8_7+Y(i,7)*A210_8_8+Y(i,8)*A210_8_9
      end do

c     Z=Z-h*[f0 F]*(B2)' 
      do i=1,m 
            Z(i,1)=Z(i,1)-h*(f0(i)*B210_1_1 + F(i,1)*gamma10)
            Z(i,2)=Z(i,2)-h*(f0(i)*B210_2_1 + F(i,2)*gamma10)
            Z(i,3)=Z(i,3)-h*(f0(i)*B210_3_1 + F(i,3)*gamma10)
            Z(i,4)=Z(i,4)-h*(f0(i)*B210_4_1 + F(i,4)*gamma10)
            Z(i,5)=Z(i,5)-h*(f0(i)*B210_5_1 + F(i,5)*gamma10)
            Z(i,6)=Z(i,6)-h*(f0(i)*B210_6_1 + F(i,6)*gamma10)
            Z(i,7)=Z(i,7)-h*(f0(i)*B210_7_1 + F(i,7)*gamma10)
            Z(i,8)=Z(i,8)-h*(f0(i)*B210_8_1 + F(i,8)*gamma10)
      end do

c     theta*Z=Z
      call sollu(m,theta,ldlu,Z(1,1),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,4),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,5),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,6),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,7),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,8),mljac,mujac,ipvt,ijob)

c
c     Y = Y-Z 
c
      do j=1,k
         do i=1,m
           Y(i,j) = Y(i,j) - Z(i,j)
         end do
      end do


      return
      end

C -------------------------------------------------------------------

      subroutine blendstep12(m,y0,f0,Y,F,h,theta,ipvt,Z,
     &                      ldlu,mljac,mujac,ijob)
c
c      Blended iteration for the 12th order method
c
      implicit none
      include 'param.dat'
      integer k
      parameter (k=10)

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision y0(m),f0(m),theta(ldlu,m),h

c     I/O parameters
      double precision Y(m,k),F(m,k)

c     Output parameters
      double precision Z(m,k)

c     Local variables
      integer i,j


C ---------------------------------------------------------------------------------------
C 12th order BIM 
C ---------------------------------------------------------------------------------------


      double precision 
     &   DA12_1_1,DA12_1_2,DA12_1_3,DA12_1_4,DA12_1_5,DA12_1_6,DA12_1_7,
     &   DA12_1_8,DA12_1_9,DA12_1_10,DA12_1_11,
     &   DA12_2_1,DA12_2_2,DA12_2_3,DA12_2_4,DA12_2_5,DA12_2_6,DA12_2_7,
     &   DA12_2_8,DA12_2_9,DA12_2_10,DA12_2_11,
     &   DA12_3_1,DA12_3_2,DA12_3_3,DA12_3_4,DA12_3_5,DA12_3_6,DA12_3_7,
     &   DA12_3_8,DA12_3_9,DA12_3_10,DA12_3_11,
     &   DA12_4_1,DA12_4_2,DA12_4_3,DA12_4_4,DA12_4_5,DA12_4_6,DA12_4_7,
     &   DA12_4_8,DA12_4_9,DA12_4_10,DA12_4_11,
     &   DA12_5_1,DA12_5_2,DA12_5_3,DA12_5_4,DA12_5_5,DA12_5_6,DA12_5_7,
     &   DA12_5_8,DA12_5_9,DA12_5_10,DA12_5_11,
     &   DA12_6_1,DA12_6_2,DA12_6_3,DA12_6_4,DA12_6_5,DA12_6_6,DA12_6_7,
     &   DA12_6_8,DA12_6_9,DA12_6_10,DA12_6_11,
     &   DA12_7_1,DA12_7_2,DA12_7_3,DA12_7_4,DA12_7_5,DA12_7_6,DA12_7_7,
     &   DA12_7_8,DA12_7_9,DA12_7_10,DA12_7_11,
     &   DA12_8_1,DA12_8_2,DA12_8_3,DA12_8_4,DA12_8_5,DA12_8_6,DA12_8_7,
     &   DA12_8_8,DA12_8_9,DA12_8_10,DA12_8_11,
     &   DA12_9_1,DA12_9_2,DA12_9_3,DA12_9_4,DA12_9_5,DA12_9_6,DA12_9_7,
     &   DA12_9_8,DA12_9_9,DA12_9_10,DA12_9_11,
     &   DA12_10_1,DA12_10_2,DA12_10_3,DA12_10_4,DA12_10_5,DA12_10_6,
     &   DA12_10_7,DA12_10_8,DA12_10_9,DA12_10_10,DA12_10_11,
     &   DB12_1_1,DB12_1_2,DB12_1_3,DB12_1_4,DB12_1_5,DB12_1_6,DB12_1_7,
     &   DB12_1_8,DB12_1_9,DB12_1_10,DB12_1_11,
     &   DB12_2_1,DB12_2_2,DB12_2_3,DB12_2_4,DB12_2_5,DB12_2_6,DB12_2_7,
     &   DB12_2_8,DB12_2_9,DB12_2_10,DB12_2_11,
     &   DB12_3_1,DB12_3_2,DB12_3_3,DB12_3_4,DB12_3_5,DB12_3_6,DB12_3_7,
     &   DB12_3_8,DB12_3_9,DB12_3_10,DB12_3_11,
     &   DB12_4_1,DB12_4_2,DB12_4_3,DB12_4_4,DB12_4_5,DB12_4_6,DB12_4_7,
     &   DB12_4_8,DB12_4_9,DB12_4_10,DB12_4_11,
     &   DB12_5_1,DB12_5_2,DB12_5_3,DB12_5_4,DB12_5_5,DB12_5_6,DB12_5_7,
     &   DB12_5_8,DB12_5_9,DB12_5_10,DB12_5_11,
     &   DB12_6_1,DB12_6_2,DB12_6_3,DB12_6_4,DB12_6_5,DB12_6_6,DB12_6_7,
     &   DB12_6_8,DB12_6_9,DB12_6_10,DB12_6_11,
     &   DB12_7_1,DB12_7_2,DB12_7_3,DB12_7_4,DB12_7_5,DB12_7_6,DB12_7_7,
     &   DB12_7_8,DB12_7_9,DB12_7_10,DB12_7_11,
     &   DB12_8_1,DB12_8_2,DB12_8_3,DB12_8_4,DB12_8_5,DB12_8_6,DB12_8_7,
     &   DB12_8_8,DB12_8_9,DB12_8_10,DB12_8_11,
     &   DB12_9_1,DB12_9_2,DB12_9_3,DB12_9_4,DB12_9_5,DB12_9_6,DB12_9_7,
     &   DB12_9_8,DB12_9_9,DB12_9_10,DB12_9_11,
     &   DB12_10_1,DB12_10_2,DB12_10_3,DB12_10_4,DB12_10_5,DB12_10_6,
     &   DB12_10_7,DB12_10_8,DB12_10_9,DB12_10_10,DB12_10_11,
     &   A212_1_1,A212_1_2,A212_1_3,A212_1_4,A212_1_5,A212_1_6,A212_1_7,
     &   A212_1_8,A212_1_9,A212_1_10,A212_1_11,
     &   A212_2_1,A212_2_2,A212_2_3,A212_2_4,A212_2_5,A212_2_6,A212_2_7,
     &   A212_2_8,A212_2_9,A212_2_10,A212_2_11,
     &   A212_3_1,A212_3_2,A212_3_3,A212_3_4,A212_3_5,A212_3_6,A212_3_7,
     &   A212_3_8,A212_3_9,A212_3_10,A212_3_11,
     &   A212_4_1,A212_4_2,A212_4_3,A212_4_4,A212_4_5,A212_4_6,A212_4_7,
     &   A212_4_8,A212_4_9,A212_4_10,A212_4_11,
     &   A212_5_1,A212_5_2,A212_5_3,A212_5_4,A212_5_5,A212_5_6,A212_5_7,
     &   A212_5_8,A212_5_9,A212_5_10,A212_5_11,
     &   A212_6_1,A212_6_2,A212_6_3,A212_6_4,A212_6_5,A212_6_6,A212_6_7,
     &   A212_6_8,A212_6_9,A212_6_10,A212_6_11,
     &   A212_7_1,A212_7_2,A212_7_3,A212_7_4,A212_7_5,A212_7_6,A212_7_7,
     &   A212_7_8,A212_7_9,A212_7_10,A212_7_11,
     &   A212_8_1,A212_8_2,A212_8_3,A212_8_4,A212_8_5,A212_8_6,A212_8_7,
     &   A212_8_8,A212_8_9,A212_8_10,A212_8_11,
     &   A212_9_1,A212_9_2,A212_9_3,A212_9_4,A212_9_5,A212_9_6,A212_9_7,
     &   A212_9_8,A212_9_9,A212_9_10,A212_9_11,
     &   A212_10_1,A212_10_2,A212_10_3,A212_10_4,A212_10_5,A212_10_6,
     &   A212_10_7,A212_10_8,A212_10_9,A212_10_10,A212_10_11,
     &   B212_1_1,B212_2_1,B212_3_1,B212_4_1,B212_5_1,B212_6_1,B212_7_1,
     &   B212_8_1,B212_9_1,B212_10_1

      parameter(
     &  DA12_1_1    =-1D0-477009404842909D0/2000000000000000D0,
     &  DA12_1_2    =1D0    +99472720252463D0/45000000000000D0,
     &  DA12_1_3    =   -104423832752463D0/20000000000000D0,
     &  DA12_1_4    =   +9994398083607D0/1250000000000D0,
     &  DA12_1_5    =   -198624609755747D0/20000000000000D0,
     &  DA12_1_6    =   +578987204267241D0/62500000000000D0,
     &  DA12_1_7    =   -63206136585249D0/10000000000000D0,
     &  DA12_1_8    =   +3823813464403D0/1250000000000D0,
     &  DA12_1_9    =   -79610832752463D0/80000000000000D0,
     &  DA12_1_10   =   +976467842623D0/5000000000000D0,
     &  DA12_1_11   =   -78691832752463D0/4500000000000000D0,
     &  DA12_2_1    =   -2813274212551D0/2929687500000D0,
     &  DA12_2_2    =   -1504054699D0/87890625000D0,
     &  DA12_2_3    =1D0    +178687703171D0/156250000000D0,
     &  DA12_2_4    =   -17253955733D0/7324218750D0,
     &  DA12_2_5    =   +137190949259D0/58593750000D0,
     &  DA12_2_6    =   -735194523679D0/366210937500D0,
     &  DA12_2_7    =   +153604208387D0/117187500000D0,
     &  DA12_2_8    =   -4522677457D0/7324218750D0,
     &  DA12_2_9    =   +23115833993D0/117187500000D0,
     &  DA12_2_10   =   -1119363137D0/29296875000D0,
     &  DA12_2_11   =   +119096296921D0/35156250000000D0,
     &  DA12_3_1    =   -4387223158569623D0/6000000000000000D0,
     &  DA12_3_2    =   -14076290383771D0/15000000000000D0,
     &  DA12_3_3    =   +45847433651313D0/20000000000000D0,
     &  DA12_3_4    =   -8091971633771D0/3750000000000D0,
     &  DA12_3_5    =   +73204095186397D0/20000000000000D0,
     &  DA12_3_6    =   -244942223059191D0/62500000000000D0,
     &  DA12_3_7    =   +84461845186397D0/30000000000000D0,
     &  DA12_3_8    =   -1752429876253D0/1250000000000D0,
     &  DA12_3_9    =   +37162883651313D0/80000000000000D0,
     &  DA12_3_10   =   -461779273473D0/5000000000000D0,
     &  DA12_3_11   =   +12525477883771D0/1500000000000000D0,
     &  DA12_4_1    =-1D0   -396187951487D0/1464843750000D0,
     &  DA12_4_2    =   +163653743573D0/175781250000D0,
     &  DA12_4_3    =   -167243587323D0/78125000000D0,
     &  DA12_4_4    =   +19779235397D0/4882812500D0,
     &  DA12_4_5    =   -105419144767D0/29296875000D0,
     &  DA12_4_6    =   +944544955011D0/244140625000D0,
     &  DA12_4_7    =   -113325741529D0/39062500000D0,
     &  DA12_4_8    =   +21327766939D0/14648437500D0,
     &  DA12_4_9    =   -9443080653D0/19531250000D0,
     &  DA12_4_10   =   +5635786799D0/58593750000D0,
     &  DA12_4_11   =   -152884212323D0/17578125000000D0,
     &  DA12_5_1    =-1D0   -5171213D0/204800000D0,
     &  DA12_5_2    =   +375871D0/4608000D0,
     &  DA12_5_3    =   -980573D0/6144000D0,
     &  DA12_5_4    =   +76277D0/384000D0,
     &  DA12_5_5    =   +456743D0/6144000D0,
     &  DA12_5_6    =1D0    +2836953D0/6400000D0,
     &  DA12_5_7    =   -2592499D0/3072000D0,
     &  DA12_5_8    =   +116713D0/384000D0,
     &  DA12_5_9    =   -2156893D0/24576000D0,
     &  DA12_5_10   =   +8271D0/512000D0,
     &  DA12_5_11   =   -640543D0/460800000D0,
     &  DA12_6_1    =   -4727199180623D0/5859375000000D0,
     &  DA12_6_2    =   -38397576209D0/58593750000D0,
     &  DA12_6_3    =   +227154597879D0/156250000000D0,
     &  DA12_6_4    =   -36961638709D0/14648437500D0,
     &  DA12_6_5    =   +123083508919D0/39062500000D0,
     &  DA12_6_6    =   -625420975389D0/244140625000D0,
     &  DA12_6_7    =1D0    +243654127213D0/117187500000D0,
     &  DA12_6_8    =   -7331573387D0/4882812500D0,
     &  DA12_6_9    =   +34452186063D0/78125000000D0,
     &  DA12_6_10   =   -4904591801D0/58593750000D0,
     &  DA12_6_11   =   +86487730543D0/11718750000000D0,
     &  DA12_7_1    =-1D0   -527365601261023D0/6000000000000000D0,
     &  DA12_7_2    =   +13219473015287D0/45000000000000D0,
     &  DA12_7_3    =   -12656585515287D0/20000000000000D0,
     &  DA12_7_4    =   +1312472696143D0/1250000000000D0,
     &  DA12_7_5    =   -72835248607009D0/60000000000000D0,
     &  DA12_7_6    =   +53134186107009D0/62500000000000D0,
     &  DA12_7_7    =   +663222376999D0/10000000000000D0,
     &  DA12_7_8    =1D0    +113677395041D0/3750000000000D0,
     &  DA12_7_9    =   -32920535515287D0/80000000000000D0,
     &  DA12_7_10   =   +906563815381D0/15000000000000D0,
     &  DA12_7_11   =   -21662785515287D0/4500000000000000D0,
     &  DA12_8_1    =-1D0   -40792014167D0/976562500000D0,
     &  DA12_8_2    =   +3448272719D0/21972656250D0,
     &  DA12_8_3    =   -1517467504D0/3662109375D0,
     &  DA12_8_4    =   +3256140146D0/3662109375D0,
     &  DA12_8_5    =   -85969177441D0/58593750000D0,
     &  DA12_8_6    =   +172929352099D0/91552734375D0,
     &  DA12_8_7    =   -14493362443D0/7324218750D0,
     &  DA12_8_8    =   +7139238634D0/3662109375D0,
     &  DA12_8_9    =   5459539781D0/39062500000D0,
     &  DA12_8_10   =   -311981201D0/2441406250D0,
     &  DA12_8_11   =   +3176506039D0/549316406250D0,
     &  DA12_9_1    =   -1883458029218509D0/2000000000000000D0,
     &  DA12_9_2    =   -1274779497593D0/5000000000000D0,
     &  DA12_9_3    =   +16125452978337D0/20000000000000D0,
     &  DA12_9_4    =   -2480966997593D0/1250000000000D0,
     &  DA12_9_5    =   +72364256949453D0/20000000000000D0,
     &  DA12_9_6    =   -308280545848359D0/62500000000000D0,
     &  DA12_9_7    =   +51140018983151D0/10000000000000D0,
     &  DA12_9_8    =   -5198771570397D0/1250000000000D0,
     &  DA12_9_9    =   +239442452978337D0/80000000000000D0,
     &  DA12_9_10   =   -979844944177D0/5000000000000D0,
     &  DA12_9_11   =   -31292283002407D0/500000000000000D0,
     &  DA12_10_1   =-1D0   -6433D0/100000D0,
     &  DA12_10_2   =   +6433D0/9000D0,
     &  DA12_10_3   =   -57897D0/16000D0,
     &  DA12_10_4   =   +2757D0/250D0,
     &  DA12_10_5   =   -45031D0/2000D0,
     &  DA12_10_6   =   +405279D0/12500D0,
     &  DA12_10_7   =   -135093D0/4000D0,
     &  DA12_10_8   =   +6433D0/250D0,
     &  DA12_10_9   =   -57897D0/4000D0,
     &  DA12_10_10  =   +6433D0/1000D0,
     &  DA12_10_11  =   -3183139D0/3600000D0)
        parameter(
     &  A212_1_1    =   477009404842909D0/2000000000000000D0,
     &  A212_1_2    =   -99472720252463D0/45000000000000D0,
     &  A212_1_3    =   104423832752463D0/20000000000000D0,
     &  A212_1_4    =   -9994398083607D0/1250000000000D0,
     &  A212_1_5    =   198624609755747D0/20000000000000D0,
     &  A212_1_6    =   -578987204267241D0/62500000000000D0,
     &  A212_1_7    =   63206136585249D0/10000000000000D0,
     &  A212_1_8    =   -3823813464403D0/1250000000000D0,
     &  A212_1_9    =   79610832752463D0/80000000000000D0,
     &  A212_1_10   =   -976467842623D0/5000000000000D0,
     &  A212_1_11   =   78691832752463D0/4500000000000000D0,
     &  A212_2_1    =   -116413287449D0/2929687500000D0,
     &  A212_2_2    =   1504054699D0/87890625000D0,
     &  A212_2_3    =   -178687703171D0/156250000000D0,
     &  A212_2_4    =   17253955733D0/7324218750D0,
     &  A212_2_5    =   -137190949259D0/58593750000D0,
     &  A212_2_6    =   735194523679D0/366210937500D0,
     &  A212_2_7    =   -153604208387D0/117187500000D0,
     &  A212_2_8    =   4522677457D0/7324218750D0,
     &  A212_2_9    =   -23115833993D0/117187500000D0,
     &  A212_2_10   =   1119363137D0/29296875000D0,
     &  A212_2_11   =   -119096296921D0/35156250000000D0,
     &  A212_3_1    =   -1612776841430377D0/6000000000000000D0,
     &  A212_3_2    =   14076290383771D0/15000000000000D0,
     &  A212_3_3    =   -45847433651313D0/20000000000000D0,
     &  A212_3_4    =   11841971633771D0/3750000000000D0,
     &  A212_3_5    =   -73204095186397D0/20000000000000D0,
     &  A212_3_6    =   244942223059191D0/62500000000000D0,
     &  A212_3_7    =   -84461845186397D0/30000000000000D0,
     &  A212_3_8    =   1752429876253D0/1250000000000D0,
     &  A212_3_9    =   -37162883651313D0/80000000000000D0,
     &  A212_3_10   =   461779273473D0/5000000000000D0,
     &  A212_3_11   =   -12525477883771D0/1500000000000000D0,
     &  A212_4_1    =   396187951487D0/1464843750000D0,
     &  A212_4_2    =   -163653743573D0/175781250000D0,
     &  A212_4_3    =   167243587323D0/78125000000D0,
     &  A212_4_4    =   -19779235397D0/4882812500D0,
     &  A212_4_5    =   134716019767D0/29296875000D0,
     &  A212_4_6    =   -944544955011D0/244140625000D0,
     &  A212_4_7    =   113325741529D0/39062500000D0,
     &  A212_4_8    =   -21327766939D0/14648437500D0,
     &  A212_4_9    =   9443080653D0/19531250000D0,
     &  A212_4_10   =   -5635786799D0/58593750000D0,
     &  A212_4_11   =   152884212323D0/17578125000000D0,
     &  A212_5_1    =   5171213D0/204800000D0,
     &  A212_5_2    =   -375871D0/4608000D0,
     &  A212_5_3    =   980573D0/6144000D0,
     &  A212_5_4    =   -76277D0/384000D0,
     &  A212_5_5    =   -456743D0/6144000D0,
     &  A212_5_6    =   -2836953D0/6400000D0,
     &  A212_5_7    =   2592499D0/3072000D0,
     &  A212_5_8    =   -116713D0/384000D0,
     &  A212_5_9    =   2156893D0/24576000D0,
     &  A212_5_10   =   -8271D0/512000D0,
     &  A212_5_11   =   640543D0/460800000D0,
     &  A212_6_1    =   -1132175819377D0/5859375000000D0,
     &  A212_6_2    =   38397576209D0/58593750000D0,
     &  A212_6_3    =   -227154597879D0/156250000000D0,
     &  A212_6_4    =   36961638709D0/14648437500D0,
     &  A212_6_5    =   -123083508919D0/39062500000D0,
     &  A212_6_6    =   625420975389D0/244140625000D0,
     &  A212_6_7    =   -243654127213D0/117187500000D0,
     &  A212_6_8    =   7331573387D0/4882812500D0,
     &  A212_6_9    =   -34452186063D0/78125000000D0,
     &  A212_6_10   =   4904591801D0/58593750000D0,
     &  A212_6_11   =   -86487730543D0/11718750000000D0,
     &  A212_7_1    =   527365601261023D0/6000000000000000D0,
     &  A212_7_2    =   -13219473015287D0/45000000000000D0,
     &  A212_7_3    =   12656585515287D0/20000000000000D0,
     &  A212_7_4    =   -1312472696143D0/1250000000000D0,
     &  A212_7_5    =   72835248607009D0/60000000000000D0,
     &  A212_7_6    =   -53134186107009D0/62500000000000D0,
     &  A212_7_7    =   -663222376999D0/10000000000000D0,
     &  A212_7_8    =   -113677395041D0/3750000000000D0,
     &  A212_7_9    =   32920535515287D0/80000000000000D0,
     &  A212_7_10   =   -906563815381D0/15000000000000D0,
     &  A212_7_11   =   21662785515287D0/4500000000000000D0,
     &  A212_8_1    =   40792014167D0/976562500000D0,
     &  A212_8_2    =   -3448272719D0/21972656250D0,
     &  A212_8_3    =   1517467504D0/3662109375D0,
     &  A212_8_4    =   -3256140146D0/3662109375D0,
     &  A212_8_5    =   85969177441D0/58593750000D0,
     &  A212_8_6    =   -172929352099D0/91552734375D0,
     &  A212_8_7    =   14493362443D0/7324218750D0,
     &  A212_8_8    =   -7139238634D0/3662109375D0,
     &  A212_8_9    =   33602960219D0/39062500000D0,
     &  A212_8_10   =   311981201D0/2441406250D0,
     &  A212_8_11   =   -3176506039D0/549316406250D0,
     &  A212_9_1    =   -116541970781491D0/2000000000000000D0,
     &  A212_9_2    =   1274779497593D0/5000000000000D0,
     &  A212_9_3    =   -16125452978337D0/20000000000000D0,
     &  A212_9_4    =   2480966997593D0/1250000000000D0,
     &  A212_9_5    =   -72364256949453D0/20000000000000D0,
     &  A212_9_6    =    308280545848359D0/62500000000000D0,
     &  A212_9_7    =   -51140018983151D0/10000000000000D0,
     &  A212_9_8    =   5198771570397D0/1250000000000D0,
     &  A212_9_9    =   -239442452978337D0/80000000000000D0,
     &  A212_9_10   =   5979844944177D0/5000000000000D0,
     &  A212_9_11   =   31292283002407D0/500000000000000D0,
     &  A212_10_1   =   6433D0/100000D0,
     &  A212_10_2   =   -6433D0/9000D0,
     &  A212_10_3   =   57897D0/16000D0,
     &  A212_10_4   =   -2757D0/250D0,
     &  A212_10_5   =   45031D0/2000D0,
     &  A212_10_6   =   -405279D0/12500D0,
     &  A212_10_7   =   135093D0/4000D0,
     &  A212_10_8   =   -6433D0/250D0,
     &  A212_10_9   =   57897D0/4000D0,
     &  A212_10_10  =   -6433D0/1000D0,
     &  A212_10_11  =   6783139D0/3600000D0)
        parameter(
     &  DB12_1_1    =   11479451563170288931D0/29536650000000000000D0,
     &  DB12_1_2    =   63885874501D0/94517280000D0,
     &  DB12_1_3    =   -57527974223D0/35286451200D0,
     &  DB12_1_4    =   2126874941D0/882161280D0,
     &  DB12_1_5    =   -5375477389D0/1960358400D0,
     &  DB12_1_6    =   5059789129D0/2205403200D0,
     &  DB12_1_7    =   -2704182763D0/1960358400D0,
     &  DB12_1_8    =   2542514071D0/4410806400D0,
     &  DB12_1_9    =   -1112711029D0/7057290240D0,
     &  DB12_1_10   =   663357529D0/26464838400D0,
     &  DB12_1_11   =   -36358589D0/21171870720D0,
     &  DB12_2_1    =   43464304639852711D0/173066308593750000D0,
     &  DB12_2_2    =   443438449D0/248107860D0,
     &  DB12_2_3    =   -3078028627D0/1722971250D0,
     &  DB12_2_4    =   313340009D0/103378275D0,
     &  DB12_2_5    =   -3589321627D0/827026200D0,
     &  DB12_2_6    =   61699387D0/13783770D0,
     &  DB12_2_7    =   -549273503D0/165405240D0,
     &  DB12_2_8    =   178915433D0/103378275D0,
     &  DB12_2_9    =   -332030663D0/551350800D0,
     &  DB12_2_10   =   156015413D0/1240539300D0,
     &  DB12_2_11   =   -59216207D0/4962157200D0,
     &  DB12_3_1    =   38675312700454505333D0/206756550000000000000D0,
     &  DB12_3_2    =   44161715831D0/26464838400D0,
     &  DB12_3_3    =   -3283281181D0/11762150400D0,
     &  DB12_3_4    =   221229581497D0/110270160000D0,
     &  DB12_3_5    =   -9249915641D0/3528645120D0, 
     &  DB12_3_6    =   328488001D0/147026880D0, 
     &  DB12_3_7    =   -24501545899D0/17643225600D0, 
     &  DB12_3_8    =   2674826191D0/4410806400D0, 
     &  DB12_3_9    =   -2067552667D0/11762150400D0, 
     &  DB12_3_10   =   160146433D0/5292967680D0, 
     &  DB12_3_11   =   -246451241D0/105859353600D0,
     &  DB12_4_1    =   148916128991770691D0/403821386718750000D0, 
     &  DB12_4_2    =   174240916D0/103378275D0, 
     &  DB12_4_3    =   -2553949D0/6891885D0, 
     &  DB12_4_4    =   112512368D0/34459425D0, 
     &  DB12_4_5    =   -4311283649D0/1531530000D0, 
     &  DB12_4_6    =   79999624D0/34459425D0, 
     &  DB12_4_7    =   -5864678D0/3828825D0, 
     &  DB12_4_8    =   4859248D0/6891885D0, 
     &  DB12_4_9    =   -7453148D0/34459425D0, 
     &  DB12_4_10   =   587116D0/14768325D0, 
     &  DB12_4_11   =   -342427D0/103378275D0,
     &  DB12_5_1    =   17710238627183D0/63515612160000D0, 
     &  DB12_5_2    =   5562130625D0/3175780608D0, 
     &  DB12_5_3    =   -935636375D0/1411458048D0, 
     &  DB12_5_4    =   2125610875D0/529296768D0, 
     &  DB12_5_5    =   -6501442375D0/2117187072D0, 
     &  DB12_5_6    =   221113406161D0/55135080000D0, 
     &  DB12_5_7    =   -6539730625D0/2117187072D0, 
     &  DB12_5_8    =   829280125D0/529296768D0, 
     &  DB12_5_9    =   -107876375D0/201636864D0, 
     &  DB12_5_10   =   349458875D0/3175780608D0, 
     &  DB12_5_11   =   -130922555D0/12703122432D0,
     &  DB12_6_1    =   84117653394000791D0/403821386718750000D0, 
     &  DB12_6_2    =   706138327D0/413513100D0, 
     &  DB12_6_3    =   -86081069D0/183783600D0, 
     &  DB12_6_4    =   24183767D0/6891885D0, 
     &  DB12_6_5    =   -615563609D0/275675400D0, 
     &  DB12_6_6    =   94460189D0/22972950D0, 
     &  DB12_6_7    =   -30823666891D0/13783770000D0, 
     &  DB12_6_8    =   32702707D0/34459425D0, 
     &  DB12_6_9    =   -1649311D0/5250960D0, 
     &  DB12_6_10   =   25613911D0/413513100D0, 
     &  DB12_6_11   =   -1836449D0/330810480D0,
     &  DB12_7_1    =   9066850209038767819D0/29536650000000000000D0, 
     &  DB12_7_2    =   1273945301D0/756138240D0, 
     &  DB12_7_3    =   -1831435457D0/5040921600D0, 
     &  DB12_7_4    =   2028629887D0/630115200D0, 
     &  DB12_7_5    =   -475125427D0/280051200D0, 
     &  DB12_7_6    =   213226433D0/63011520D0, 
     &  DB12_7_7    =   -25983713D0/56010240D0, 
     &  DB12_7_8    =   8396816521D0/15752880000D0, 
     &  DB12_7_9    =   -1310408951D0/5040921600D0, 
     &  DB12_7_10   =   175300267D0/3780691200D0, 
     &  DB12_7_11   =   -58184383D0/15122764800D0,
     &  DB12_8_1    =   49560837081155911D0/173066308593750000D0, 
     &  DB12_8_2    =   76433504D0/44304975D0, 
     &  DB12_8_3    =   -18894032D0/34459425D0, 
     &  DB12_8_4    =   385150592D0/103378275D0, 
     &  DB12_8_5    =   -54124624D0/20675655D0, 
     &  DB12_8_6    =   6293440D0/1378377D0, 
     &  DB12_8_7    =   -165253856D0/103378275D0, 
     &  DB12_8_8    =   250726016D0/103378275D0, 
     &  DB12_8_9    =   -8968348841D0/13783770000D0, 
     &  DB12_8_10   =   3948064D0/62026965D0, 
     &  DB12_8_11   =   -1974032D0/310134825D0,
     &  DB12_9_1    =   2538784308052589777D0/9845550000000000000D0, 
     &  DB12_9_2    =   14976693257D0/8821612800D0, 
     &  DB12_9_3    =   -326374859D0/784143360D0, 
     &  DB12_9_4    =   4913986843D0/1470268800D0, 
     &  DB12_9_5    =   -10987159111D0/5881075200D0, 
     &  DB12_9_6    =   862239689D0/245044800D0, 
     &  DB12_9_7    =   -2973275233D0/5881075200D0, 
     &  DB12_9_8    =   443339993D0/294053760D0, 
     &  DB12_9_9    =   4141950047D0/3920716800D0, 
     &  DB12_9_10   =   -52868905381D0/220540320000D0, 
     &  DB12_9_11   =   -411440411D0/35286451200D0,
     &  DB12_10_1   =   80335D0/299376D0, 
     &  DB12_10_2   =   132875D0/74844D0, 
     &  DB12_10_3   =   -80875D0/99792D0, 
     &  DB12_10_4   =   28375D0/6237D0, 
     &  DB12_10_5   =   -24125D0/5544D0, 
     &  DB12_10_6   =   89035D0/12474D0, 
     &  DB12_10_7   =   -24125D0/5544D0, 
     &  DB12_10_8   =   28375D0/6237D0, 
     &  DB12_10_9   =   -80875D0/99792D0, 
     &  DB12_10_10  =   132875D0/74844D0, 
     &  DB12_10_11  =   -8769811D0/23388750D0,
     &  B212_1_1    =   -5169648083607D0/50000000000000D0,
     &  B212_2_1    =   97697971D0/6103515625D0,
     &  B212_3_1    =   4558075961257D0/50000000000000D0,
     &  B212_4_1    =   -8992156761D0/97656250000D0,
     &  B212_5_1    =   -45031D0/5120000D0,
     &  B212_6_1    =   1614856691D0/24414062500D0,
     &  B212_7_1    =   -1513503946143D0/50000000000000D0,
     &  B212_8_1    =   -1333129889D0/97656250000D0,
     &  B212_9_1    =   872716997593D0/50000000000000D0,
     &  B212_10_1   =   0D0)
C ---------------------------------------------------------------------------------------

c     Z=[y0 Y]*(DA)'
      do i=1,m
         Z(i,1)=y0(i) *DA12_1_1  +Y(i,1) *DA12_1_2  +Y(i,2)*DA12_1_3+
     &            Y(i,3)*DA12_1_4  +Y(i,4) *DA12_1_5  +Y(i,5)*DA12_1_6+
     &            Y(i,6)*DA12_1_7  +Y(i,7) *DA12_1_8  +Y(i,8)*DA12_1_9+
     &            Y(i,9)*DA12_1_10 +Y(i,10)*DA12_1_11
         Z(i,2)=y0(i) *DA12_2_1  +Y(i,1) *DA12_2_2  +Y(i,2)*DA12_2_3+ 
     &            Y(i,3)*DA12_2_4  +Y(i,4) *DA12_2_5  +Y(i,5)*DA12_2_6+
     &            Y(i,6)*DA12_2_7  +Y(i,7) *DA12_2_8  +Y(i,8)*DA12_2_9+
     &            Y(i,9)*DA12_2_10 +Y(i,10)*DA12_2_11
           Z(i,3)=y0(i) *DA12_3_1  +Y(i,1) *DA12_3_2  +Y(i,2)*DA12_3_3+
     &            Y(i,3)*DA12_3_4  +Y(i,4) *DA12_3_5  +Y(i,5)*DA12_3_6+
     &            Y(i,6)*DA12_3_7  +Y(i,7) *DA12_3_8  +Y(i,8)*DA12_3_9+
     &            Y(i,9)*DA12_3_10 +Y(i,10)*DA12_3_11
           Z(i,4)=y0(i) *DA12_4_1  +Y(i,1) *DA12_4_2  +Y(i,2)*DA12_4_3+
     &            Y(i,3)*DA12_4_4  +Y(i,4) *DA12_4_5  +Y(i,5)*DA12_4_6+
     &            Y(i,6)*DA12_4_7  +Y(i,7) *DA12_4_8  +Y(i,8)*DA12_4_9+
     &            Y(i,9)*DA12_4_10 +Y(i,10)*DA12_4_11
           Z(i,5)=y0(i) *DA12_5_1  +Y(i,1) *DA12_5_2  +Y(i,2)*DA12_5_3+
     &            Y(i,3)*DA12_5_4  +Y(i,4) *DA12_5_5  +Y(i,5)*DA12_5_6+
     &            Y(i,6)*DA12_5_7  +Y(i,7) *DA12_5_8  +Y(i,8)*DA12_5_9+
     &            Y(i,9)*DA12_5_10 +Y(i,10)*DA12_5_11
           Z(i,6)=y0(i) *DA12_6_1  +Y(i,1) *DA12_6_2  +Y(i,2)*DA12_6_3+
     &            Y(i,3)*DA12_6_4  +Y(i,4) *DA12_6_5  +Y(i,5)*DA12_6_6+
     &            Y(i,6)*DA12_6_7  +Y(i,7) *DA12_6_8  +Y(i,8)*DA12_6_9+
     &            Y(i,9)*DA12_6_10 +Y(i,10)*DA12_6_11
           Z(i,7)=y0(i) *DA12_7_1  +Y(i,1) *DA12_7_2  +Y(i,2)*DA12_7_3+
     &            Y(i,3)*DA12_7_4  +Y(i,4) *DA12_7_5  +Y(i,5)*DA12_7_6+
     &            Y(i,6)*DA12_7_7  +Y(i,7) *DA12_7_8  +Y(i,8)*DA12_7_9+
     &            Y(i,9)*DA12_7_10 +Y(i,10)*DA12_7_11
           Z(i,8)=y0(i) *DA12_8_1  +Y(i,1) *DA12_8_2  +Y(i,2)*DA12_8_3+
     &            Y(i,3)*DA12_8_4  +Y(i,4) *DA12_8_5  +Y(i,5)*DA12_8_6+
     &            Y(i,6)*DA12_8_7  +Y(i,7) *DA12_8_8  +Y(i,8)*DA12_8_9+
     &            Y(i,9)*DA12_8_10 +Y(i,10)*DA12_8_11
           Z(i,9)=y0(i) *DA12_9_1  +Y(i,1) *DA12_9_2  +Y(i,2)*DA12_9_3+
     &            Y(i,3)*DA12_9_4  +Y(i,4) *DA12_9_5  +Y(i,5)*DA12_9_6+
     &            Y(i,6)*DA12_9_7  +Y(i,7) *DA12_9_8  +Y(i,8)*DA12_9_9+
     &            Y(i,9)*DA12_9_10 +Y(i,10)*DA12_9_11
          Z(i,10)=y0(i) *DA12_10_1 +Y(i,1) *DA12_10_2 +Y(i,2)*DA12_10_3+
     &            Y(i,3)*DA12_10_4 +Y(i,4) *DA12_10_5 +Y(i,5)*DA12_10_6+
     &            Y(i,6)*DA12_10_7 +Y(i,7) *DA12_10_8 +Y(i,8)*DA12_10_9+
     &            Y(i,9)*DA12_10_10+Y(i,10)*DA12_10_11
      end do

c     Z=Z-h*[f0 F]*(dB)'
      do i=1,m 
         Z(i,1) = Z(i,1)-
     &         h*(f0(i) *DB12_1_1  +F(i,1) *DB12_1_2  +F(i,2)*DB12_1_3+
     &            F(i,3)*DB12_1_4  +F(i,4) *DB12_1_5  +F(i,5)*DB12_1_6+
     &            F(i,6)*DB12_1_7  +F(i,7) *DB12_1_8  +F(i,8)*DB12_1_9+
     &            F(i,9)*DB12_1_10 +F(i,10)*DB12_1_11)
         Z(i,2) = Z(i,2)-
     &         h*(f0(i) *DB12_2_1  +F(i,1) *DB12_2_2  +F(i,2)*DB12_2_3+
     &            F(i,3)*DB12_2_4  +F(i,4) *DB12_2_5  +F(i,5)*DB12_2_6+
     &            F(i,6)*DB12_2_7  +F(i,7) *DB12_2_8  +F(i,8)*DB12_2_9+
     &            F(i,9)*DB12_2_10 +F(i,10)*DB12_2_11)
         Z(i,3) = Z(i,3)-
     &         h*(f0(i) *DB12_3_1  +F(i,1) *DB12_3_2  +F(i,2)*DB12_3_3+
     &            F(i,3)*DB12_3_4  +F(i,4) *DB12_3_5  +F(i,5)*DB12_3_6+
     &            F(i,6)*DB12_3_7  +F(i,7) *DB12_3_8  +F(i,8)*DB12_3_9+
     &            F(i,9)*DB12_3_10 +F(i,10)*DB12_3_11)
         Z(i,4) = Z(i,4)-
     &         h*(f0(i) *DB12_4_1  +F(i,1) *DB12_4_2  +F(i,2)*DB12_4_3+
     &            F(i,3)*DB12_4_4  +F(i,4) *DB12_4_5  +F(i,5)*DB12_4_6+
     &            F(i,6)*DB12_4_7  +F(i,7) *DB12_4_8  +F(i,8)*DB12_4_9+
     &            F(i,9)*DB12_4_10 +F(i,10)*DB12_4_11)
         Z(i,5) = Z(i,5)-
     &         h*(f0(i) *DB12_5_1  +F(i,1) *DB12_5_2  +F(i,2)*DB12_5_3+
     &            F(i,3)*DB12_5_4  +F(i,4) *DB12_5_5  +F(i,5)*DB12_5_6+
     &            F(i,6)*DB12_5_7  +F(i,7) *DB12_5_8  +F(i,8)*DB12_5_9+
     &            F(i,9)*DB12_5_10 +F(i,10)*DB12_5_11)
         Z(i,6) = Z(i,6)-
     &         h*(f0(i) *DB12_6_1  +F(i,1) *DB12_6_2  +F(i,2)*DB12_6_3+
     &            F(i,3)*DB12_6_4  +F(i,4) *DB12_6_5  +F(i,5)*DB12_6_6+
     &            F(i,6)*DB12_6_7  +F(i,7) *DB12_6_8  +F(i,8)*DB12_6_9+
     &            F(i,9)*DB12_6_10 +F(i,10)*DB12_6_11)
         Z(i,7) = Z(i,7)-
     &         h*(f0(i) *DB12_7_1  +F(i,1) *DB12_7_2  +F(i,2)*DB12_7_3+
     &            F(i,3)*DB12_7_4  +F(i,4) *DB12_7_5  +F(i,5)*DB12_7_6+
     &            F(i,6)*DB12_7_7  +F(i,7) *DB12_7_8  +F(i,8)*DB12_7_9+
     &            F(i,9)*DB12_7_10 +F(i,10)*DB12_7_11)
         Z(i,8) = Z(i,8)-
     &         h*(f0(i) *DB12_8_1  +F(i,1) *DB12_8_2  +F(i,2)*DB12_8_3+
     &            F(i,3)*DB12_8_4  +F(i,4) *DB12_8_5  +F(i,5)*DB12_8_6+
     &            F(i,6)*DB12_8_7  +F(i,7) *DB12_8_8  +F(i,8)*DB12_8_9+
     &            F(i,9)*DB12_8_10 +F(i,10)*DB12_8_11)
         Z(i,9) = Z(i,9)-
     &         h*(f0(i) *DB12_9_1  +F(i,1) *DB12_9_2  +F(i,2)*DB12_9_3+
     &            F(i,3)*DB12_9_4  +F(i,4) *DB12_9_5  +F(i,5)*DB12_9_6+
     &            F(i,6)*DB12_9_7  +F(i,7) *DB12_9_8  +F(i,8)*DB12_9_9+
     &            F(i,9)*DB12_9_10 +F(i,10)*DB12_9_11)
        Z(i,10) = Z(i,10)-
     &         h*(f0(i) *DB12_10_1 +F(i,1) *DB12_10_2 +F(i,2)*DB12_10_3+
     &            F(i,3)*DB12_10_4 +F(i,4) *DB12_10_5 +F(i,5)*DB12_10_6+
     &            F(i,6)*DB12_10_7 +F(i,7) *DB12_10_8 +F(i,8)*DB12_10_9+
     &            F(i,9)*DB12_10_10+F(i,10)*DB12_10_11)
      end do

c     theta*Z=Z
      call sollu(m,theta,ldlu,Z(1,1),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,4),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,5),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,6),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,7),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,8),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,9),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,10),mljac,mujac,ipvt,ijob)

c     Z= Z+[y0 Y]*(A2)'
      do i=1,m
       Z(i,1)=Z(i,1)+
     &            y0(i)* A212_1_1  +Y(i,1) *A212_1_2  +Y(i,2)*A212_1_3+
     &           Y(i,3)* A212_1_4  +Y(i,4) *A212_1_5  +Y(i,5)*A212_1_6+
     &           Y(i,6)* A212_1_7  +Y(i,7) *A212_1_8  +Y(i,8)*A212_1_9+
     &           Y(i,9)* A212_1_10 +Y(i,10)*A212_1_11
         Z(i,2)=Z(i,2)+
     &            y0(i)* A212_2_1  +Y(i,1) *A212_2_2  +Y(i,2)*A212_2_3+ 
     &           Y(i,3)* A212_2_4  +Y(i,4) *A212_2_5  +Y(i,5)*A212_2_6+
     &           Y(i,6)* A212_2_7  +Y(i,7) *A212_2_8  +Y(i,8)*A212_2_9+
     &           Y(i,9)* A212_2_10 +Y(i,10)*A212_2_11
         Z(i,3)=Z(i,3)+
     &            y0(i)* A212_3_1  +Y(i,1) *A212_3_2  +Y(i,2)*A212_3_3+
     &           Y(i,3)* A212_3_4  +Y(i,4) *A212_3_5  +Y(i,5)*A212_3_6+
     &           Y(i,6)* A212_3_7  +Y(i,7) *A212_3_8  +Y(i,8)*A212_3_9+
     &           Y(i,9)* A212_3_10 +Y(i,10)*A212_3_11
         Z(i,4)=Z(i,4)+
     &            y0(i)* A212_4_1  +Y(i,1) *A212_4_2  +Y(i,2)*A212_4_3+
     &           Y(i,3)* A212_4_4  +Y(i,4) *A212_4_5  +Y(i,5)*A212_4_6+
     &           Y(i,6)* A212_4_7  +Y(i,7) *A212_4_8  +Y(i,8)*A212_4_9+
     &           Y(i,9)* A212_4_10 +Y(i,10)*A212_4_11
         Z(i,5)=Z(i,5)+
     &            y0(i)* A212_5_1  +Y(i,1) *A212_5_2  +Y(i,2)*A212_5_3+
     &           Y(i,3)* A212_5_4  +Y(i,4) *A212_5_5  +Y(i,5)*A212_5_6+
     &           Y(i,6)* A212_5_7  +Y(i,7) *A212_5_8  +Y(i,8)*A212_5_9+
     &           Y(i,9)* A212_5_10 +Y(i,10)*A212_5_11
         Z(i,6)=Z(i,6)+
     &            y0(i)* A212_6_1  +Y(i,1) *A212_6_2  +Y(i,2)*A212_6_3+
     &           Y(i,3)* A212_6_4  +Y(i,4) *A212_6_5  +Y(i,5)*A212_6_6+
     &           Y(i,6)* A212_6_7  +Y(i,7) *A212_6_8  +Y(i,8)*A212_6_9+
     &           Y(i,9)* A212_6_10 +Y(i,10)*A212_6_11
         Z(i,7)=Z(i,7)+
     &            y0(i)* A212_7_1  +Y(i,1) *A212_7_2  +Y(i,2)*A212_7_3+
     &           Y(i,3)* A212_7_4  +Y(i,4) *A212_7_5  +Y(i,5)*A212_7_6+
     &           Y(i,6)* A212_7_7  +Y(i,7) *A212_7_8  +Y(i,8)*A212_7_9+
     &           Y(i,9)* A212_7_10 +Y(i,10)*A212_7_11
         Z(i,8)=Z(i,8)+
     &            y0(i)* A212_8_1  +Y(i,1) *A212_8_2  +Y(i,2)*A212_8_3+
     &           Y(i,3)* A212_8_4  +Y(i,4) *A212_8_5  +Y(i,5)*A212_8_6+
     &           Y(i,6)* A212_8_7  +Y(i,7) *A212_8_8  +Y(i,8)*A212_8_9+
     &           Y(i,9)* A212_8_10 +Y(i,10)*A212_8_11
       Z(i,9)=Z(i,9)+
     &            y0(i)* A212_9_1  +Y(i,1) *A212_9_2  +Y(i,2)*A212_9_3+
     &           Y(i,3)* A212_9_4  +Y(i,4) *A212_9_5  +Y(i,5)*A212_9_6+
     &           Y(i,6)* A212_9_7  +Y(i,7) *A212_9_8  +Y(i,8)*A212_9_9+
     &           Y(i,9)* A212_9_10 +Y(i,10)*A212_9_11
      Z(i,10)=Z(i,10)+
     &            y0(i)* A212_10_1 +Y(i,1) *A212_10_2 +Y(i,2)*A212_10_3+
     &           Y(i,3)* A212_10_4 +Y(i,4) *A212_10_5 +Y(i,5)*A212_10_6+
     &           Y(i,6)* A212_10_7 +Y(i,7) *A212_10_8 +Y(i,8)*A212_10_9+
     &           Y(i,9)* A212_10_10+Y(i,10)*A212_10_11
      end do

c     Z=Z-h*[f0 F]*(B2)'
      do i=1,m 
            Z(i,1)=Z(i,1) -h*(f0(i)*B212_1_1  + F(i,1)*gamma12)
            Z(i,2)=Z(i,2) -h*(f0(i)*B212_2_1  + F(i,2)*gamma12)
            Z(i,3)=Z(i,3) -h*(f0(i)*B212_3_1  + F(i,3)*gamma12)
            Z(i,4)=Z(i,4) -h*(f0(i)*B212_4_1  + F(i,4)*gamma12)
            Z(i,5)=Z(i,5) -h*(f0(i)*B212_5_1  + F(i,5)*gamma12)
            Z(i,6)=Z(i,6) -h*(f0(i)*B212_6_1  + F(i,6)*gamma12)
            Z(i,7)=Z(i,7) -h*(f0(i)*B212_7_1  + F(i,7)*gamma12)
            Z(i,8)=Z(i,8) -h*(f0(i)*B212_8_1  + F(i,8)*gamma12)
            Z(i,9)=Z(i,9) -h*(f0(i)*B212_9_1  + F(i,9)*gamma12)
           Z(i,10)=Z(i,10)-h*(f0(i)*B212_10_1 + F(i,10)*gamma12)
      end do

c     theta*Z=Z
      call sollu(m,theta,ldlu,Z(1,1),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,4),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,5),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,6),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,7),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,8),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,9),mljac,mujac,ipvt,ijob)
      call sollu(m,theta,ldlu,Z(1,10),mljac,mujac,ipvt,ijob)

c
c     Y = Y-Z 
c
      do j=1,k
         do i=1,m
           Y(i,j) = Y(i,j) - Z(i,j)
         end do
      end do


      return
      end

C -------------------------------------------------------------------

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 
C     LOCALERR, ERRDOWN AND ERRUP
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      subroutine localerr4(m,f0,F,h,Z,scal,nerr,nerrup,nlinsys,
     &                     theta,ipvt,ldlu,mljac,mujac,ijob)
c
c     Local error estimate for the method of order 4 
c
      implicit none
      include 'param.dat'

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision f0(m),F(m,3),theta(ldlu,m),h,scal(m)


c     Output parameters
      double precision Z(m,3),nerr,nerrup

c     I/O parameters
      integer nlinsys

c     Local variables
      integer i
      double precision VMAX2

      do i=1,m
           Z(i,1)=h*(PSI4_1*f0(i) +PSI4_2*F(i,1)-
     &               PSI4_2*F(i,2)-PSI4_1*F(i,3))
           Z(i,2)=Z(i,1)
      end do
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      
      do i=1,m
            Z(i,3)=Z(i,1)-Z(i,2)
      end do
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      VMAX2 = gamma4*VMAX4_2
      do i=1,m
           Z(i,2)=VMAX4_1*Z(i,2)
           Z(i,3)=VMAX2*Z(i,3)
      end do
    
      call norm(m,2,scal,Z(1,2),nerr,nerrup)

      nlinsys = nlinsys + 2

      return
      end

C -------------------------------------------------------------------

      subroutine localerr6(m,f0,F,h,Z,scal,nerr,nerrup,nlinsys,
     &                     theta,ipvt,ldlu,mljac,mujac,ijob)
c
c     Local error estimate for the method of order 6 
c
      implicit none
      include 'param.dat'

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision f0(m),F(m,4),theta(ldlu,m),h,scal(m)

c     Output parameters
      double precision Z(m,3),nerr,nerrup

c     I/O parameters
      integer nlinsys

c     Local variables
      integer i
      double precision VMAX2

      do i=1,m
           Z(i,1)=h*(PSI6_1*f0(i) +PSI6_2*F(i,1)+
     &               PSI6_3*F(i,2)+PSI6_2*F(i,3)+
     &               PSI6_1*F(i,4))
           Z(i,2)=Z(i,1)
      end do
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)

      do i=1,m
            Z(i,3)=2d0*Z(i,1)-Z(i,2)
      end do
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      do i=1,m
            Z(i,3)=Z(i,1)-Z(i,3)
      end do
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      VMAX2 = gamma6*VMAX6_2
      do i=1,m
           Z(i,2)=VMAX6_1*Z(i,2)
           Z(i,3)=VMAX2*Z(i,3)
      end do
      call norm(m,2,scal,Z(1,2),nerr,nerrup)
      nlinsys = nlinsys + 3
      return
      end

C -------------------------------------------------------------------
      
      subroutine errdown6(m,f0,F,h,Z,scal,nerrdown,nlinsys,
     &                    qinf,theta,ipvt,ldlu,mljac,mujac,ijob)
c
c     Estimate of the error for the lower-order method,  k=4 
c
      implicit none
      include 'param.dat'

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision f0(m),F(m,3),theta(ldlu,m),h,scal(m)
      logical qinf

c     Output parameters
      double precision Z(m),nerrdown

c     I/O parameters
      integer nlinsys

c     Local variables
      integer i
      double precision VMAX2
      if (qinf) then 
         VMAX2 = gamma6*VMAX4_2
      else
         VMAX2 = VMAX4_1
      end if
      do i=1,m
           Z(i)=VMAX2*h*(PSI4_1*f0(i) +PSI4_2*F(i,1)-
     &                   PSI4_2*F(i,2)-PSI4_1*F(i,3))
      end do
      call sollu(m,theta,ldlu,Z,mljac,mujac,ipvt,ijob)
      call norm(m,1,scal,Z,nerrdown,VMAX2)

      nlinsys = nlinsys + 1

      return
      end

C -------------------------------------------------------------------
C -------------------------------------------------------------------

      subroutine localerr8(m,f0,F,h,Z,scal,nerr,nerrup,nlinsys,
     &                     theta,ipvt,ldlu,mljac,mujac,ijob)
c
c     Local error estimate for the method of order 8
c
      implicit none
      include 'param.dat'

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision f0(m),F(m,6),theta(ldlu,m),h,scal(m)

c     Output parameters
      double precision Z(m,3),nerr,nerrup

c     I/O parameters
      integer nlinsys


c     Local variables
      integer i
      double precision VMAX2

      do i=1,m
           Z(i,1)=h*(PSI8_1*f0(i) +PSI8_2*F(i,1)+
     &               PSI8_3*F(i,2)+PSI8_4*F(i,3)+
     &               PSI8_3*F(i,4)+PSI8_2*F(i,5)+
     &               PSI8_1*F(i,6))     
           Z(i,2)=Z(i,1)
      end do
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)

      do i=1,m
            Z(i,3)=2d0*Z(i,1)-Z(i,2)
      end do
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      do i=1,m
            Z(i,3)=Z(i,1)-Z(i,3)
      end do
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      VMAX2 = gamma8*VMAX8_2
      do i=1,m
           Z(i,2)=VMAX8_1*Z(i,2)
           Z(i,3)=VMAX2*Z(i,3)
      end do
      call norm(m,2,scal,Z(1,2),nerr,nerrup)
      nlinsys = nlinsys + 3
      return
      end

C -------------------------------------------------------------------
      
      subroutine errdown8(m,f0,F,h,Z,scal,nerrdown,nlinsys,
     &                    qinf,theta,ipvt,ldlu,mljac,mujac,ijob)
c
c     Estimate of the error for the lower-order method,  k=6
c
      implicit none
      include 'param.dat'

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision f0(m),F(m,4),theta(ldlu,m),h,scal(m)
      logical qinf

c     Output parameters
      double precision Z(m),nerrdown

c     I/O parameters
      integer nlinsys

c     Local variables
      integer i
      double precision VMAX2
      if (qinf) then 
        VMAX2 = gamma8*VMAX6_2
      else
        VMAX2 = VMAX6_1
      end if 
      do i=1,m
           Z(i)=VMAX2*h*(PSI6_1*f0(i) +PSI6_2*F(i,1)+
     &                   PSI6_3*F(i,2)+PSI6_2*F(i,3)+
     &                   PSI6_1*F(i,4))
      end do
      call sollu(m,theta,ldlu,Z,mljac,mujac,ipvt,ijob)
      call norm(m,1,scal,Z,nerrdown,VMAX2)

      nlinsys = nlinsys + 1

      return
      end

C -------------------------------------------------------------------
C -------------------------------------------------------------------

      subroutine localerr10(m,f0,F,h,Z,scal,nerr,nerrup,
     &                     nlinsys,theta,ipvt,ldlu,mljac,mujac,ijob)
c
c     Local error estimate for the method of order 10
c
      implicit none
      include 'param.dat'

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision f0(m),F(m,8),theta(ldlu,m),h,scal(m)

c     Output parameters
      double precision Z(m,3),nerr,nerrup

c     I/O parameters
      integer nlinsys

c     Local variables
      integer i
      double precision VMAX2

      do i=1,m
           Z(i,1)=h*(PSI10_1*f0(i) +PSI10_2*F(i,1)+
     &               PSI10_3*F(i,2)+PSI10_4*F(i,3)+
     &               PSI10_5*F(i,4)+PSI10_4*F(i,5)+
     &               PSI10_3*F(i,6)+PSI10_2*F(i,7)+
     &               PSI10_1*F(i,8))
           Z(i,2)=Z(i,1)
      end do
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)

      do i=1,m
            Z(i,3)=2d0*Z(i,1)-Z(i,2)
      end do
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      do i=1,m
            Z(i,3)=Z(i,1)-Z(i,3)
      end do
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)

      VMAX2 = gamma10*VMAX10_2
      do i=1,m
           Z(i,2)=VMAX10_1*Z(i,2)
           Z(i,3)=VMAX2*Z(i,3)
      end do
      call norm(m,2,scal,Z(1,2),nerr,nerrup)
      nlinsys = nlinsys + 3
      return
      end

C -------------------------------------------------------------------
      
      subroutine errdown10(m,f0,F,h,Z,scal,nerrdown,nlinsys,
     &                    qinf,theta,ipvt,ldlu,mljac,mujac,ijob)
c
c     Estimate of the error for the lower-order method,  k=8
c
      implicit none
      include 'param.dat'

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision f0(m),F(m,6),theta(ldlu,m),h,scal(m)
      logical qinf

c     Output parameters
      double precision Z(m),nerrdown

c     I/O parameters
      integer nlinsys

c     Local variables
      integer i
      double precision VMAX2
      if (qinf) then 
        VMAX2 = gamma10*VMAX8_2
      else
        VMAX2 = VMAX8_1
      end if
      do i=1,m
           Z(i)=VMAX2*h*(PSI8_1*f0(i) +PSI8_2*F(i,1)+
     &                   PSI8_3*F(i,2)+PSI8_4*F(i,3)+
     &                   PSI8_3*F(i,4)+PSI8_2*F(i,5)+
     &                   PSI8_1*F(i,6))     
      end do
      call sollu(m,theta,ldlu,Z,mljac,mujac,ipvt,ijob)
      call norm(m,1,scal,Z,nerrdown,VMAX2)

      nlinsys = nlinsys + 1

      return
      end
 
C -------------------------------------------------------------------
C -------------------------------------------------------------------

      subroutine localerr12(m,f0,F,h,Z,scal,nerr,nerrup,
     &                     nlinsys,theta,ipvt,ldlu,mljac,mujac,ijob)
c
c     Local error estimate for the method of order 12
c
      implicit none
      include 'param.dat'

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision f0(m),F(m,10),theta(ldlu,m),h,scal(m)

c     Output parameters
      double precision Z(m,3),nerr,nerrup

c     I/O parameters
      integer nlinsys

c     Local variables
      integer i
      double precision VMAX2
      do i=1,m
           Z(i,1)=h*(PSI12_1*f0(i) +PSI12_2*F(i,1)+
     &               PSI12_3*F(i,2)+PSI12_4*F(i,3)+
     &               PSI12_5*F(i,4)+PSI12_6*F(i,5)+
     &               PSI12_5*F(i,6)+PSI12_4*F(i,7)+
     &               PSI12_3*F(i,8)+PSI12_2*F(i,9)+
     &               PSI12_1*F(i,10))
           Z(i,2)=Z(i,1)
      end do
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)

      do i=1,m
            Z(i,3)=2d0*Z(i,1)-Z(i,2)
      end do
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)
      do i=1,m
            Z(i,3)=Z(i,1)-Z(i,3)
      end do
      call sollu(m,theta,ldlu,Z(1,3),mljac,mujac,ipvt,ijob)

      VMAX2 = gamma12*VMAX12_2
      do i=1,m
           Z(i,2)=VMAX12_1*Z(i,2)
           Z(i,3)=VMAX2*Z(i,3)
      end do
      call norm(m,2,scal,Z(1,2),nerr,nerrup)
      nlinsys = nlinsys + 3
      return
      end

C -------------------------------------------------------------------
      
      subroutine errdown12(m,f0,F,h,Z,scal,nerrdown,nlinsys,
     &                    qinf,theta,ipvt,ldlu,mljac,mujac,ijob)
c
c     Estimate of the error for the lower-order method,  k=10
c
      implicit none
      include 'param.dat'

c     Input parameters
      integer m,ipvt(m),ldlu,mljac,mujac,ijob
      double precision f0(m),F(m,8),theta(ldlu,m),h,scal(m)
      logical qinf

c     Output parameters
      double precision Z(m),nerrdown

c     I/O parameters
      integer nlinsys

c     Local variables
      integer i
      double precision VMAX2
      if(qinf) then
         VMAX2 = gamma12*VMAX10_2
      else
         VMAX2 = VMAX10_1
      end if
      do i=1,m
           Z(i)=VMAX2*h*(PSI10_1*f0(i) +PSI10_2*F(i,1)+
     &                   PSI10_3*F(i,2)+PSI10_4*F(i,3)+
     &                   PSI10_5*F(i,4)+PSI10_4*F(i,5)+
     &                   PSI10_3*F(i,6)+PSI10_2*F(i,7)+
     &                   PSI10_1*F(i,8))
      end do
      call sollu(m,theta,ldlu,Z,mljac,mujac,ipvt,ijob)
      call norm(m,1,scal,Z,nerrdown,VMAX2)

      nlinsys = nlinsys + 1

      return
      end

C -------------------------------------------------------------------
C ---------------------- SUBROUTINE ERRUP ------------------------------

      subroutine errup(m,k,ord_ind,Z,h,h0,h00,vup,nerrup,scal,
     &                 theta,ipvt,ldlu,mljac,mujac,ijob)
c
c     Estimate of the higher-order method, in case of order reduction
c     
      implicit none 
c     Input parameter
      integer m,k,ipvt(m),ldlu,mljac,mujac,ijob,ord_ind
      double precision Z(m,k+2),theta(ldlu,m),scal(m),
     &                 h,h0,h00,vup

c     Output parameter 
      double precision nerrup

c     Local variables
      integer i
      double precision dbk,dh,dh0,cp,al
      

      goto (10,20,20,20) ord_ind

10    continue
      dh = h/h0
      cp = vup/dble(k)
      dbk= dble(k+1)
      dh = dh**dbk
      do i=1,m
         Z(i,2)=cp * ( Z(i,1)-Z(i,k+1)*dh)
      end do
      goto 30

20    continue
      dbk = dble(k+1)
      al  = h/h0
      dh  = al**dbk
      dh0 = (h/h00)**dbk
      cp  = 2d0*vup/dble(k*k)*h/(h+h0)
      do i=1,m
        Z(i,2) = cp*(   (   Z(i,1)  - dh *Z(i,k+1)) - 
     &               al*(dh*Z(i,k+1)- dh0*Z(i,k+2)) )
      end do
      goto 30

30    continue
      call sollu(m,theta,ldlu,Z(1,2),mljac,mujac,ipvt,ijob)
      call norm(m,1,scal,Z(1,2),nerrup,dh)

      return

      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 
C     EXTRAPOLA, DIFFDIV AND NORM
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE EXTRAPOLA (M,K,KNEW,H0,H,YNEW,DD)
C
C     Initial profile by extrapolation. 
C

      IMPLICIT NONE

C     Input parameters
      INTEGER M,K,KNEW
      DOUBLE PRECISION H0,H,DD(K+1,M)

C     Output parameters
      DOUBLE PRECISION YNEW(M,KNEW)

C     Local variables
      INTEGER I,J,L
      DOUBLE PRECISION DT,RATH

      RATH = (H/H0)

      DO I=1,M
C     Evaluation of the interpolating polynomial at the points of the new block        
        DO L=1,KNEW
           DT  = dble(L)*RATH
           YNEW(I,L)=DD(K+1,I)
           DO J=K,1,-1
              DT = DT+1d0
              YNEW(I,L)=YNEW(I,L)*DT +DD(J,I)
           END DO
        END DO
      END DO

      RETURN
      END

      SUBROUTINE DIFFDIV(M,K,Y0,Y,DD)
C     Compute the divided differences of the interpolating polynomial

      IMPLICIT NONE

C     Input parameters
      INTEGER M,K
	DOUBLE PRECISION Y0(M),Y(M,K)

C     Output parameters
      DOUBLE PRECISION DD(K+1,M)
	
	      
C     Local variables
      INTEGER I,J,L
	DOUBLE PRECISION DT

      DO I=1,M
          DD(1,I)=Y0(I)
          DO J=1,K
            DD(J+1,I)=Y(I,J)
          END DO
          DO J=2,K+1
            DT = dble(J-1)
            DO L=K+1,J,-1
               DD(L,I)=(DD(L,I)-DD(L-1,I))/DT
            END DO
          END DO
      END DO
      RETURN
	END

      SUBROUTINE NORM(M,K,SCAL,ERR,NERR,NERRUP)
C
C     Used norm
C
      IMPLICIT NONE 

C     INPUT PARAMETER
      INTEGER M,K
      DOUBLE PRECISION SCAL(M),ERR(M,K)

C     OUTPUT PARAMETER
      DOUBLE PRECISION NERR,NERRUP

C     LOCAL VARIABLES
      DOUBLE PRECISION NERR0
      INTEGER I,J

      NERR = 0d0
      DO J=1,K-1
         NERR0=0d0
         DO I=1,M
           NERR0 = NERR0 + ( ERR(I,J)*SCAL(I))**2
         END DO
         NERR = DMAX1(NERR,NERR0)
      END DO

      NERRUP=0d0
      DO I=1,M
        NERRUP = NERRUP + ( ERR(I,K)*SCAL(I))**2
      END DO
      NERR = DMAX1(NERR,NERRUP)

      NERR   = DSQRT(NERR/DBLE(M))
      NERRUP = DSQRT(NERRUP/DBLE(M))
      RETURN
      END   

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 
C     LINEAR ALGEBRA
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 
C The routines DEC, DECB, SOL and SOLB have been replaced by interface routines
C that call LINPACK routines instead.
C
C     SUBROUTINE DEC
C
      SUBROUTINE DEC (N, NDIM, A, IP, IER)
      implicit none
C VERSION REAL DOUBLE PRECISION
      INTEGER N,NDIM,IP,IER
      DOUBLE PRECISION A
      DIMENSION A(NDIM,N), IP(N)
C-----------------------------------------------------------------------
C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
C  INPUT..
C     N = ORDER OF MATRIX.
C     NDIM = DECLARED DIMENSION OF ARRAY  A .
C     A = MATRIX TO BE TRIANGULARIZED.
C  OUTPUT..
C     A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
C     A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
C     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
C     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
C     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
C           SINGULAR AT STAGE K.
C  USE  SOL  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
C  DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
C  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
C
C  REFERENCE..
C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER,
C     C.A.C.M. 15 (1972), P. 274.
C-----------------------------------------------------------------------
c     
c     The following replace the original DEC routine with a call to the
c     the LAPACK routine DGETRF
      integer M
      M = N
      call DGETRF( M, N, A, ndim, IP, ier )
      if (ier < 0) then
         write(6,*) 'Error in DGETRF call via DEC, INFO = ', ier
         stop
      end if
      RETURN
C----------------------- END OF SUBROUTINE DEC -------------------------
      END
C
C     SUBROUTINE SOL
C
      SUBROUTINE SOL (N, NDIM, A, B, IP)
C VERSION REAL DOUBLE PRECISION
      implicit none
      INTEGER N,NDIM,IP
      DOUBLE PRECISION A,B
      DIMENSION A(NDIM,N), B(N), IP(N)
C-----------------------------------------------------------------------
C  SOLUTION OF LINEAR SYSTEM, A*X = B .
C  INPUT..
C    N = ORDER OF MATRIX.
C    NDIM = DECLARED DIMENSION OF ARRAY  A .
C    A = TRIANGULARIZED MATRIX OBTAINED FROM DEC.
C    B = RIGHT HAND SIDE VECTOR.
C    IP = PIVOT VECTOR OBTAINED FROM DEC.
C  DO NOT USE IF DEC HAS SET IER .NE. 0.
C  OUTPUT..
C    B = SOLUTION VECTOR, X .
C-----------------------------------------------------------------------
C
C  The original code has been replaced with a call to the LAPACK routine
C  DGETRS
      INTEGER            INFO
      call DGETRS( 'N', N, 1, A, ndim, IP, B, N, INFO )
      if (info < 0) then
         write(6,*) 'Error in DGETRS call via DGESL, INFO = ', info
      stop
      end if
      RETURN
C----------------------- END OF SUBROUTINE SOL -------------------------
      END


      subroutine sollu(n,a,lda,b,ml,mu,ipvt,ijob)
      integer lda,n,ipvt(n),ijob
      double precision a(lda,n),b(n)

      goto (1,2) ijob

1     call sol(n,n,a,b,ipvt)
      return

2     call solb(n,lda,a,ml,mu,b,ipvt)
      return
      end

      subroutine declu(n,a,lda,ml,mu,ipvt,ijob,info)
      integer lda,n,ipvt(n),info,ijob
      double precision a(lda,n)

      goto(1,2) ijob

   1  call dec(n,n,a,ipvt,info)
      return
    
   2  call decb(n,lda,a,ml,mu,ipvt,info)
      return

      end

    
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccca
c
C     SUBROUTINE DECB
C
      SUBROUTINE DECB (N, NDIM, A, ML, MU, IP, IER)
      implicit none
      REAL*8 A
      integer n, ndim, ml, mu, ip, ier
      DIMENSION A(NDIM,N), IP(N)
C-----------------------------------------------------------------------
C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED
C  MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU
C  INPUT..
C     N       ORDER OF THE ORIGINAL MATRIX A.
C     NDIM    DECLARED DIMENSION OF ARRAY  A.
C     A       CONTAINS THE MATRIX IN BAND STORAGE.   THE COLUMNS
C                OF THE MATRIX ARE STORED IN THE COLUMNS OF  A  AND
C                THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
C                ML+1 THROUGH 2*ML+MU+1 OF  A.
C     ML      LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
C     MU      UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
C  OUTPUT..
C     A       AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
C                THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
C     IP      INDEX VECTOR OF PIVOT INDICES.
C     IP(N)   (-1)**(NUMBER OF INTERCHANGES) OR O .
C     IER     = 0 IF MATRIX A IS NONSINGULAR, OR  = K IF FOUND TO BE
C                SINGULAR AT STAGE K.
C  USE  SOLB  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
C  DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N)  WITH MD=ML+MU+1.
C  IF IP(N)=O, A IS SINGULAR, SOLB WILL DIVIDE BY ZERO.
C
C  REFERENCE..
C     THIS IS A MODIFICATION OF
C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER,
C     C.A.C.M. 15 (1972), P. 274.
C-----------------------------------------------------------------------
C  This routine is replaced by a call to the LAPACK routine DGBTRF
      call DGBTRF( N, N, ML, MU, a, ndim, IP, ier )  
      if (ier < 0) then
         write(6,*) 'Error in DGBTRF call via DECB, INFO = ', ier
         stop
      end if
      RETURN
C----------------------- END OF SUBROUTINE DECB ------------------------
      END
C
C     SUBROUTINE SOLB
C
      SUBROUTINE SOLB (N, NDIM, A, ML, MU, B, IP)
      implicit none
      REAL*8 A,B
      integer n, ndim, ml, mu, ip
      DIMENSION A(NDIM,N), B(N), IP(N)
C-----------------------------------------------------------------------
C  SOLUTION OF LINEAR SYSTEM, A*X = B .
C  INPUT..
C    N      ORDER OF MATRIX A.
C    NDIM   DECLARED DIMENSION OF ARRAY  A .
C    A      TRIANGULARIZED MATRIX OBTAINED FROM DECB.
C    ML     LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
C    MU     UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
C    B      RIGHT HAND SIDE VECTOR.
C    IP     PIVOT VECTOR OBTAINED FROM DECB.
C  DO NOT USE IF DECB HAS SET IER .NE. 0.
C  OUTPUT..
C    B      SOLUTION VECTOR, X .
C-----------------------------------------------------------------------
C  This routine is replaced by the LAPARC routine DGBTRS 
      integer info
      call DGBTRS( 'N', N, ML, MU, 1, a, ndim, IP, B, n, info )
      if (info < 0) then
         write(6,*) 'Error in DGBTRS call via DGBSL, INFO = ', info
         stop
      end if
      RETURN
C----------------------- END OF SUBROUTINE SOLB ------------------------
      END
