c===================== Edition from 27.04.2007 =========================
c
            subroutine RiSolve(F_Y          [Reference],
     *                         Param_Model  [Reference],
     *                         F_Risolv     [Reference],
     *                         I_Risolv     [Reference],
     *                         H_Risolv     [Reference],
     *                         YY_Risolv    [Reference],
     *                         U1_Risolv    [Reference],
     *                         U2_Risolv    [Reference],
     *                         LIST         [Reference],
     *                         YYstep       [Reference])
c 
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      include 'max_dimension.fi' ! NEQ - число уравнений
                
       real*8    F_Y(0:NEQ+1)     ! начальное приближение стартового решения
c                                ! из файла
	
	integer*4 I_Risolv(0:20)
 	real*8    F_Risolv(0:20)
 	
	real*8    H_Risolv (0:NEQ+1)! H - абсол.значения эталонных шагов 
	                            !по параметру 
	real*8    YY_Risolv(0:NEQ+1)! YY - масштабы для переменных 
	
	real*8    U1_Risolv(0:NEQ+1)  ! U1    массивы U1,U2 - область  
	real*8    U2_Risolv(0:NEQ+1)  ! U2    определения решения

	real*8    YYstep(0:NEQ)  
	
	include 'Param_common.fi'   ! параметры модели

	real*8 Param_Model(1:Nparam)
	common /W_Ypositive/ I_Ypositive ! контроль положительности решения
      REAL*8 U(1:NEQ+1)   
c
	real*8 Ystep(1:NEQ+4)
c
      COMMON /ADD_R/    S(Nparam)
      COMMON /Par_Mod/  S0(Nparam)
c
        common /WI_Max/ I_end, I_Max       !!!!!!!!!!!!!!!!
        common /WF_Max/ F_Max              !!!!!!!!!!!!!!!!
c
c динамические массивы
	real*8 YBIF, YQS
	allocatable YBIF(:,:) ! YBIF(I_Risolv(12),NEQ+2)
        allocatable YQS (:,:) ! YQS (I_Risolv(12),NEQ+2)
c
      common /bifur_turn_steps/ Nbif, Nturn, Nsys
c
      include 'RiSolve_DLLEXPORT.fi'
!DEC$ IF DEFINED (make_DLL)
      include 'RiSolve_dll.fi'
!DEC$ ENDIF      
c      
c-----------------------------------------------------------------------
c
      I_END=0 
c
      if (LIST.eq.0) then
      open(unit=33,file='ResTab_old.bin',status='unknown',
     #     BUFFERED = 'YES',BLOCKSIZE=20480,FORM='BINARY')
      ifail=0
      end if
c
c     Basic data input 
c
        I_Ypositive=0
	Nalpha=I_Risolv(1)    ! Nalpha - the parameter for continuation method
	NDIM=NEQ
        Nsys=NEQ
	NPAR=Nalpha
c
c Model parameters for common blockc
c
          p(:)=Param_Model(:)
          S(1:Nparam)=p(1:Nparam)
	    S0(1:Nparam)=p(1:Nparam)
c
c     Constants
c
      QMIN= F_Risolv(1)     ! QMIN
      QMAX= F_Risolv(2)     ! QMAX
      HSTEP=F_Risolv(3)     ! HSTEP
      QSOLV=F_Risolv(4)     ! qsolv = alpha
      H00=  F_Risolv(5)     ! Hpar

      U(1:NDIM)=F_Y(1:NDIM) ! x(i) initial approximation of start solution
 
      NP1=NDIM+1
      U(NP1)=S(NPAR) 
c
c     Standard data input 
      ACCEL=   F_Risolv(6)    ! 
      EXACT=   F_Risolv(7)    ! 
      EMIN=    F_Risolv(8)    ! 
      EMAX=    F_Risolv(9)    ! 
      STAB=    F_Risolv(10)   ! 
      EXACT0=  F_Risolv(11)   ! 
      IPARAM  =I_Risolv(2)    ! 
      Jmatr   =I_Risolv(3)    ! 
      LIMIT   =I_Risolv(4)    ! 
      LIMIT0  =I_Risolv(5)    ! 
      IPREDICT=I_Risolv(6)    ! 
      NDEX    =I_Risolv(7)    ! 
      Nbif    =I_Risolv(11)   !
      Nturn   =I_Risolv(12)   !
      Nstep   =I_Risolv(13)   !
      Jacob   =I_Risolv(15)   !   
      IHTM    =I_Risolv(16)   !    
      INTER   =I_Risolv(17)   !    
      H_Risolv(NP1)=H00

      NEND=Nstep

c задание массивов
	allocate ( YBIF(Nturn,NEQ+2) )
	allocate ( YQS (Nturn,NEQ+2) )
c
c*********************************************************************
        Jacob=0 ! The numerical differentiation is absent !!!

c      open (unit=32, mode='write',file='error0.DAT',status='unknown',
c     *      BUFFERED = 'YES',BLOCKSIZE=20480)
c	rewind(32)
c      write(32,*) ' '
c         write(32,*) (U(I),I=1,NP1) !!!!!!!!!
c      close(32)  

	CALL SNEQP(NDIM,NDEX,IPARAM,Jmatr,Jacob,NEND,limit,limit0,
     *           ipredict,hstep,accel,exact,emin,emax,
     *           NPAR,LIST,KBIF,ISOLV,ifail,INTER,QSOLV,QMIN,QMAX,
     *           STAB,EXACT0,U(1:NP1),
     *           H_Risolv (1:NP1), YY_Risolv(1:NP1),
     *           U1_Risolv(1:NEQ+1), U2_Risolv(1:NEQ+1),
     *           Ystep,YQS,YBIF)
      YYstep(0)=Ystep(NEQ+1)
	YYstep(1:NEQ)=Ystep(1:NEQ)
c
         if (ifail.gt.0) I_Risolv(20)=ifail
	
	   if (LIST.eq.1) then
         open(unit=33,file='ResTab_old.bin',status='REPLACE',
     #        BUFFERED = 'YES',BLOCKSIZE=20480,FORM='BINARY')
	   else
	   open(unit=33,file='ResTab_old.bin',status='OLD',
     #        access='append', BUFFERED = 'YES',BLOCKSIZE=20480,
     #        FORM='BINARY')
	   end if
c
      write(33) Ystep(:) !
      close(33)
c
      CALL error('Risolve')
c
      InPut_I_Risolv_20=I_Risolv(20)            
      I_Risolv(20)=ifail
      
c--------- Таблицы ---------------------------------------------------      
         LIMSTEP=LIST
	   if (InPut_I_Risolv_20.ge.0) then
            I_TABST=1
	   else
	    I_TABST=0 
	   end if

         CALL TABST (NDIM,KBIF,ISOLV,KHBP,LIMSTEP,NPAR,IHTM,
     *               STAB,QSOLV,Ystep,YQS,YBIF,I_TABST)

c     deallocate ( YBIF ) !!!!!!!!!!!!!!!!!
c     deallocate ( YQS )  !!!!!!!!!!!!!!!!!
      I_Risolv(7)= NDEX       ! NDEX=NEQ+1
c Ошибки
      CALL error('Risolve')
c
      return
      END
c--------------------------------------------------------------------------------------
c     Tables for results
c
      SUBROUTINE TABST(NDIM,KBIF,KSOL,KHBP,LIST,NPAR,IHTM,
     *                 STAB,QSOLV,Ystep,YQS,YBIF,I_TABST)
      implicit integer*4 (i-n)
      implicit real*8 (A-H,O-Z)
	  include 'number_equations.fi'
      include 'Param_common.fi'
      common /bifur_turn_steps/ Nbif, Nturn, Nsys 

      REAL*8  YQS (Nturn,NDIM+2)
      REAL*8  YBIF(Nturn,NDIM+2)
      REAL*8  Ystep(NDIM+4)
      REAL*8  HB(Nbif*2)
      REAL*8  ZB(Nbif*2)
      DIMENSION NB(Nbif*2)

      CHARACTER*8     chu, Defx(NEQ)
      CHARACTER*14    chq*14
      CHARACTER*17    RM0(NEQ+1)
      CHARACTER*35    DATLEN,TIMLEN
      COMMON /DATA/   DATLEN,TIMLEN     
C----------------- OUTPUT the table -------------------
      NP1=NDIM+1
         if (LIST.eq.1) then
  	DO k=1,NDIM
	if(k.lt.10) then
	write(Defx(k),'(a,I1,a)') '  x',k,'    '
	else
      if(k.lt.100) then
	write(Defx(k),'(a,I2,a)') '  x',k,'   '
	else
      if(k.lt.1000) then
	write(Defx(k),'(a,I3,a)') '  x',k,'  '
	else
	write(Defx(k),'(a,I4,a)') '  x',k,' '
	end if
	end if
	end if
	END DO 
         open(unit=120,file='results.dat',status='REPLACE',
     #        BUFFERED = 'YES',BLOCKSIZE=20480)
	   write(120,998) 'alfa',(Defx(k),k=1,ndim)
	   else
	   open(unit=120,file='results.dat',status='OLD',access='append',
     #        BUFFERED = 'YES',BLOCKSIZE=20480)
	   end if
	write(120,999) Ystep(np1),(Ystep(k),k=1,ndim)
        close(unit=120)  

      if (I_TABST.eq.1)  then 

c       BASIC TABLES !

        LIMSTEP=LIST
        CALL bas_tables (NDIM,KHBP,LIMSTEP,NPAR,IHTM,NB,
     #                   STAB,Ystep,HB,ZB)
c ----------------------------------------------------------
      OPEN(UNIT=1,FILE='Bifurc_mult.dat',BUFFERED='YES',
     #     BLOCKSIZE=20480)
c ----------------------------------------------------------
      WRITE(1,906) DATLEN
      WRITE(1,906) TIMLEN
      if(KBIF.gt.0) then
	write(1,905) '  Table 2. Specific points of "turning" type on st.
     &s.d.'
       DO I=1,KBIF
	 Ns=NINT(YBIF(I,NP1+1))
       WRITE(1,'(/,4x,a,I4)') 'Number of turning point:',I
       WRITE(1,'(4x,a,e15.8)') 'Parameter value:',YBIF(I,NP1)
       WRITE(1,'(4x,a,I4)') 'Number of step at the turning point:',Ns
       write(1,'(4x,a)') 'Component values:'
       write(1,22) (YBIF(I,J),J=1,NDIM)
       write(1,115) ('-',k=1,70)
       END DO
      end if
	close(1)

       OPEN(UNIT=1,FILE='Bifurc_stabl.dat',BUFFERED='YES',
     #     BLOCKSIZE=20480)
       WRITE(1,906) DATLEN
       WRITE(1,906) TIMLEN

       if((STAB.gt.0.d0).and.(KHBP.gt.0)) then

	write(1,905) '  Table 3. Points of stability change on st.s.d.'

	DO L=1,KHBP
       WRITE(1,'(/,4x,a,I4)') 'Number of point of stability change:',L
       WRITE(1,'(4x,a,e15.8)') 'Parameter value::',HB(L)
	 if(ZB(L).lt.0.d0) then
	 WRITE(1,'(4x,a)') '(with change of determinant sign)'
	 else
	 WRITE(1,'(4x,a)') '(without change of determinant sign)'
	 end if
       WRITE(1,'(4x,a,I4)') 'Number of step at the point of stability 
     #change:',NB(L)
       write(1,115) ('-',k=1,70)
	END DO
      end if
      close(1)

      OPEN(UNIT=1,FILE='Multiplicity.dat',BUFFERED='YES',
     #     BLOCKSIZE=20480)
      WRITE(1,906) DATLEN
      WRITE(1,906) TIMLEN

      if(KSOL.gt.0) then
c
        if(STAB.gt.0.d0) then
        call reqsolv(NDIM,KSOL,NPAR,STAB,QSOLV,YQS)
        else
        write(1,'(/,1x,a,I6)') 'The number of solutions =',KSOL
	write(1,905) 'Table 4. Multiplicity solutions on st.s.d. '
	write(1,'(5x,a,I4,a,e10.4,/)') 'at the given value of model parame
     &ter p(Npar)=Qsolv, Npar=',NPAR,', Qsolv=',qsolv

        DO J=1,KSOL
          WRITE(1,'(/,4x,a,I4)') 'Number of solution:',J
          write(1,'(4x,a)') 'Component values:'
          write(1,22) (YQS(J,I),I=1,NDIM)
          write(1,115) ('-',k=1,70)
        END DO
        end if
      else
        write(1,'(/,1x,a)') 'The solutions are absent at Q=QSOLV !'
      end if
      CLOSE(UNIT=1)
c ----------------------------------------------------------
      OPEN(UNIT=1,FILE='Tab_Exl.mlt',BUFFERED='YES',
     #     BLOCKSIZE=20480)
        write(1,'(I6)') KSOL
      if(KSOL.gt.0) then
        DO J=1,KSOL
         write(1,'(I6)') J
         DO I=1,NDIM
          write(1,'(e16.8)') YQS(J,I)
         END DO
        END DO
      end if 
      CLOSE(UNIT=1)
c ----------------------------------------------------------
      close(32)   
      end if !!!!!!!!!!!!!!!!!!!!!!!!!!! (I_TABST=1)

22    FORMAT(2X,100e16.8)
115   FORMAT(70a)
904   FORMAT(1x,a,/)
905   FORMAT(/,10x,a,/)
906   FORMAT(A)
995   format(6x,a,7x,1000(4x,a8,5x))
998   format(6x,a,6x,1000(4x,a8,4x))
999   format(1000(SP,e16.8))
1000  format(I4,1000(SP,e16.8))
1001  format(a17,1000(a17))

c ERRORS !
      CALL error('TABST')
c           
      RETURN
      END
c*******************************************************************
      subroutine SNEQP(N,NDEX,IPARAM,Jmatr,Jacob,NEND,LIMIT,LIMIT0,
     *                 IPREDICT,HSTEP,ACCEL,EXACT,EMIN,EMAX,
     *                 NPAR,LIST,KBIF,NSOLV,ifail,INTER,QSOLV,QMIN,
     *                 QMAX,STAB,EXACT0,U,H,YY,U1,U2,Ystep,YQS,YBIF)
      implicit real*8 (A-H,O-Z)
      implicit integer*4 (i-n)
c
	include 'number_equations.fi'
	include 'param_model.fi'

      common /bifur_turn_steps/ Nbif, Nturn, Nsys
      real*8 YQS (Nturn,N+2)
      REAL*8 YBIF(Nturn,N+2)
c
      character*8 chu, DATE, TIME
      CHARACTER*35    DATLEN,TIMLEN
       
      REAL*8 Ystep(NEQ+4) 
      REAL*8 U00(NEQ+1),
     *       F  (NEQ+1),
     *       Y1 (NEQ+1),
     *       Y2 (NEQ+1),
     *       W1 (NEQ+1),
     *       W2 (NEQ+1),
     *       H  (N+1),
     *       YY (N+1), 
     *       U  (N+1), 
     *       U1 (N+1),
     *       U2 (N+1)

	DIMENSION IA(Nturn) 
      COMMON /ADD_R/ S(Nparam)
      COMMON /DATA/   DATLEN,TIMLEN     

      DATA DATE/' DATE : '/
      DATA TIME/' TIME : '/

      allocatable FM(:,:)
c
      if (LIST.gt.NEND) ifail=5
      if (ifail.gt.0) return 
c
c Восстановление с предыдущего шага
	iWork=0
	if (LIST.gt.0) then
	NReturn=1
	call One_Step(NReturn,Nsys,Nturn,
     *        N,NDEX,IPARAM,Jmatr,NEND,LIMIT,LIMIT0,ifail,
     *        IPREDICT,NPAR,LIST,KBIF,NSOLV,
     *        ISTART, IA(1:Nturn),NP1,N0,
     *        HSTEP,ACCEL,EXACT,EMIN,EMAX,H1,
     *        QSOLV,QMIN,QMAX,STAB,QAPPA,EXACT0,
     *        H  (1:NEQ+1),
     *        YY (1:NEQ+1), 
     *        U(1:NEQ+1), 
     *        U1(1:NEQ+1),
     *        U2(1:NEQ+1),
     *        Ystep(1:NEQ+4), YQS(1:Nturn,1:NEQ+2),
     *        YBIF(1:Nturn,1:NEQ+2),
     *        F  (1:NEQ+1),
     *        Y1 (1:NEQ+1),
     *        Y2 (1:NEQ+1),
     *        W1 (1:NEQ+1),
     *        W2 (1:NEQ+1))
 	goto 200
	end if
c
c---------  Program constants  ------------------------------------
      NP1=N+1
      N0 =N
c
      NSOLV=0
      LIST=0
      KBIF=0
      ifail=-1

	CALL GETDAT(JY,JM,JD)
        WRITE(DATLEN,'(1X,A,I2.2,1H-,I2.2,1H-,I4.4)') DATE,JD,JM,JY
        call gettim(ih,im,is,is0)
        WRITE(TIMLEN,'(1X,A,I2,A,I2,A,I2)') TIME,ih,':',im,':',is
c-----------------------------------------------------------------------
c Файл runtime error     
      open (unit=32, mode='write', file='error.DAT', status='unknown',
     *      BUFFERED = 'YES',BLOCKSIZE=20480)
	rewind(32)
      write(32,*) ' '
      write(32,*) ' version 27.04.07  '
      write(32,*) ' '
      WRITE(32,906) DATLEN
      WRITE(32,906) TIMLEN

        DO I=1,N
         if(U(I).lt.U1(I).or.U(I).gt.U2(I)) ifail=3
        END DO
	if(S(NPAR).lt.QMIN.or.S(NPAR).gt.QMAX) ifail=4
        DO I=1,NP1
         if(H(I).lt.0.d0) ifail=12
        END DO
	if(ifail.gt.0) return
c
c     Newton's method
c
      U00(1:NP1)=U(1:NP1)

      IF (Jmatr.eq.1) THEN
      allocate ( FM(Nparam,NEQ)	)
      ndex00=ndex
      DO J=1,Nparam

       Y1(1:N)=U00(1:N)
       Y1(NP1)=S(J)
       ifail0=-1
	 ndex0=ndex00
      if(J.ne.NPAR) then
       ISTAB=0
      else
       ISTAB=1
      end if 
      CALL NEWTON0 (N,J,ndex0,limit0,jacob,ifail0,exact,
     #              Y1,Y2,U1,U2,DET,STAB,QAPPA0,ISTAB)
	DO k=1,N
	 FM(J,k)=dabs(Y2(k))
	END DO
	if(J.eq.NPAR) then
	U(1:NP1)=Y1(1:NP1)
        F(1:NP1)=Y2(1:NP1)
	 ifail=ifail0
	  ndex=ndex0
          DETB=DET
          QAPPA=QAPPA0
	 end if
	END DO
         call sensFM(N,NPAR,FM)
         deallocate(FM)
	ELSE
         U(NP1)=S(NPAR)
         ISTAB=1
         CALL NEWTON0 (N,NPAR,ndex,limit0,jacob,ifail,exact,
     #                 U,F,U1,U2,DETB,STAB,QAPPA,ISTAB)
        END IF
        if (ifail.gt.0) then
         LIST=1
         return
        end if
      GOTO 200
1000  CONTINUE
      CALL NEWTON(N,NPAR,LIST,ndex,limit,jacob,ipredict,
     *            hstep,H1,exact,emin,ifail,
     *            H,U1,U2,U,Y1,Y2,W1,W2,F,DETB,STAB,QAPPA) 

       W0=W2(NP1)
       F0=F(NP1)
       Y0=Y2(NP1)
       U0=U(NP1)
       LL=DSIGN(1d0,F0)-DSIGN(1d0,W0)
       if(LL.NE.0) 
     #       CALL BIFUR (N,U,F,Y2,W2,NDEX,LIST,KBIF,YBIF)

        R1=Y0-QSOLV
        R2=U0-QSOLV
        LL=DSIGN(1d0,R2)-DSIGN(1d0,R1)
        if(LL.NE.0) then
         CALL SOLV_Q(N,U,F,Y2,W2,NDEX,NSOLV,QSOLV,YQS)
         IA(NSOLV)=NDEX
       end if
c
c      Q in  [QMIN,QMAX] ?
c
       IF(U0.GT.QMAX.OR.U0.LT.QMIN) then
        if(U0.GT.QMAX) then
         QLIMIT=QMAX
        else
         QLIMIT=QMIN
        end if
        HH=U0-Y0
        TT=(QLIMIT-Y0)/HH
        T1=1d0-TT
        T2=T1*T1
        T0=TT*TT
        T12=(1.d0+2.d0*TT)*T2
        T32=(3.d0-2.d0*TT)*T0
        HHW=HH*TT*T2/W0 
        HHF=HH*T0*T1/F0
        U(1:N)=Y2(1:N)*T12+U(1:N)*T32+HHW*W2(1:N)-HHF*F(1:N)
        U(NP1)=QLIMIT
        ifail=0
       END IF

200   continue
	iWork=iWork+1
      if (iWork.lt.2) LIST=LIST+1 
c+++++++++++++++++
      if (iWork.eq.2) then
     
      NReturn=0

	call One_Step(NReturn,Nsys,Nturn,
     *        N,NDEX,IPARAM,Jmatr,NEND,LIMIT,LIMIT0,ifail,
     *        IPREDICT,NPAR,LIST,KBIF,NSOLV,
     *        ISTART, IA(1:Nturn),NP1,N0,
     *        HSTEP,ACCEL,EXACT,EMIN,EMAX,H1,
     *        QSOLV,QMIN,QMAX,STAB,QAPPA,EXACT0,
     *        H  (1:NEQ+1),
     *        YY (1:NEQ+1), 
     *        U(1:NEQ+1), 
     *        U1(1:NEQ+1),
     *        U2(1:NEQ+1),
     *        Ystep(1:NEQ+4), YQS(1:Nturn,1:NEQ+2),
     *        YBIF(1:Nturn,1:NEQ+2),
     *        F  (1:NEQ+1),
     *        Y1 (1:NEQ+1),
     *        Y2 (1:NEQ+1),
     *        W1 (1:NEQ+1),
     *        W2 (1:NEQ+1))
	return
	end if
c++++++++++++++++++++
	DO I=1,NP1
	 Ystep(I)=U(I)
	END DO
	Ystep(N+2)=NDEX
        Ystep(N+3)=DETB
        Ystep(N+4)=QAPPA-STAB
c
        if (LIST.eq.NEND)  ifail=5
        if(KBIF.eq.Nturn)  ifail=9
        if(NSOLV.eq.Nturn) ifail=10
        if(INTER.gt.0)     ifail=11
        if (ifail.ge.0) goto 100
c
      ISTART=NDEX
c
      IF(LIST.EQ.1) THEN

c     Choice of new current parameter on start

      if (IPARAM.lt.2) then
       NN=NP1-IPARAM
       HMAX=0.d0
        DO I=1,NN
         IF(DABS(F(I)).GT.HMAX) THEN
          HMAX=DABS(F(I))
          NDEX=I
         END IF
        END DO

       FN=F(NDEX)
       F(1:NP1)=F(1:NP1)/FN
      END IF

c    Calculating the start step and initial approximation 

      DO I=1,NP1
       if(DABS(HSTEP*F(I)).gt.H(I)) HSTEP=H(I)/DABS(F(I)) 
      END DO
      SGNH=DSIGN(1d0,HSTEP)
      DIRECT=DSIGN(1d0,HSTEP*F(NP1))
      if(SGNH*DIRECT.lt.0d0) HSTEP=-HSTEP      
   
      Y2(1:NP1)=U(1:NP1)
      W2(1:NP1)=F(1:NP1)
      U(1:NP1)=Y2(1:NP1)+HSTEP*W2(1:NP1)
      ELSE
c
c     Choice of new current parameter
c
      IF(IPARAM.LT.2) THEN
c
       NN=NP1-IPARAM
       HMAX=0.d0
        DO I=1,NN
         FMAX=DABS(U(I)-Y2(I)/YY(I))
         IF(FMAX.GT.HMAX) THEN
          HMAX=FMAX
          NDEX=I
         END IF
        END DO
 
       FN=F(NDEX)
       F(1:NP1)=F(1:NP1)/FN
      END IF
c
c    Calculation of step and initial approximation 
c
c     Set of the current step
c
      H1=U(NDEX)-Y2(NDEX)
        if(DABS(H1).lt.EMIN*H(NDEX)) then
          ifail=7
          goto 100
        end if
         HSTEP=ACCEL*H1
         DO I=1,NP1
          HH=EMAX*H(I)
          if(DABS(HSTEP*F(I)).gt.HH) HSTEP=HH/DABS(F(I))
         END DO
        SGNH=DSIGN(1d0,U(ISTART)-Y2(ISTART))
        DIRECT=DSIGN(1.d0,HSTEP*F(ISTART))
        if(SGNH*DIRECT.lt.0d0) HSTEP=-HSTEP 

      WN=W2(NDEX)
      Y1(1:NP1)=Y2(1:NP1)
      W1(1:NP1)=W2(1:NP1)/WN
      Y2(1:NP1)=U(1:NP1)
      W2(1:NP1)=F(1:NP1)
c
c    Calculation of initial approximation 
c
         CALL APPROX (N,IPREDICT,HSTEP,H1,U,F,Y1,Y2,W1,W2)
         END IF

      GOTO 1000
 100  CONTINUE
      if(ifail.ne.-1) then
        call gettim(ih1,im1,is1,is01)
        itime=(ih1-ih)*60+im1-im
        write(32,'(a,I3,a)') 
     #  ' Time of calculation=',itime,' (min.)'
      write(32,'(a,I6)') 'The number of all solutions=',NSOLV
      NSOLV0=NSOLV
      if (NSOLV.gt.0) then 
       ifail0=ifail
         DO J=1,NSOLV
         U(1:N)=YQS(J,1:N)            
   	  U(NP1)=QSOLV
	  ndex0=IA(J)
          ifail=-1
          ISTAB=1
         CALL NEWTON0 (N,NPAR,ndex0,limit0,jacob,ifail,
     #                 exact,U,F,U1,U2,DETB,STAB,QAPPA,ISTAB)
            
         YQS(J,1:N)=U(1:N)            
	   YQS(J,NP1)=QSOLV
           if(STAB.gt.0d0) then
            YQS(J,N+2)=QAPPA-STAB
           else
            YQS(J,N+2)=0d0
           end if
         END DO
       ifail=ifail0
      end if
c
      if(NSOLV.gt.1) then
      J=0
10    J=J+1
      U(1:N)=YQS(J,1:N)
      K=J        
20    K=K+1
      SUM=0d0
      DO I=1,N
       SUM=SUM+(YQS(K,I)-U(I))**2
      END DO 
      Qerr=DSQRT(SUM)
      if(Qerr.lt.EXACT0) then      
       if(K.lt.NSOLV) then
        DO L=K,NSOLV-1
         DO I=1,N+2 
          YQS(L,I)=YQS(L+1,I)
         END DO 
        END DO
       end if
       NSOLV=NSOLV-1
       K=K-1
      end if
      if(K.lt.NSOLV) goto 20
      if(J.lt.NSOLV) goto 10
      end if
      end if
c
c Ошибки
      CALL error('SNEQP')
906   FORMAT(A)
c         
      RETURN
      END
c-------------------------------------------------------------
      SUBROUTINE GAUSS(N,B,F,IGAUSS,LU_dec)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (i-n)
	include 'number_equations.fi'
      REAL*8 P(NEQ)
      REAL*8 B(N,N),F(N+1)
      DIMENSION LG(NEQ), MG(NEQ)
c
      IGAUSS=0
      EPS=1.d-305
      N1=N-1
c
         IF (LU_dec.EQ.0) GOTO 25
c
      DO K=1,N1
      C0=0.d+00
         DO I=K,N
            DO J=K,N
               IF (DABS(B(I,J)).LT.DABS(C0)) cycle
               C0=B(I,J)
            L0=I
            LG(K)=I
            M0=J
            MG(K)=J
            END DO
         END DO
         IF (DABS(C0).LT.EPS) THEN
          IGAUSS=1
          GOTO 100
         END IF 


         P(1:N)=B(1:N,M0)
         B(1:N,M0)=B(1:N,K)
         B(1:N,K)=P(1:N)
         P(K:N)=B(L0,K:N)
         B(L0,K:N)=B(K,K:N)
         B(K,K)=C0
         KK=K+1
         DO J=KK,N
         B(K,J)=P(J)/C0
             DO I=KK,N
             B(I,J)=B(I,J)-B(I,K)*B(K,J)
             END DO
         END DO
      END DO
c
         IF(DABS(B(N,N)).LT.EPS) THEN
          IGAUSS=1
          goto 100
         END IF

   25 CONTINUE
         DO K=1,N1
         LK=LG(K)
         F0=F(LK)
         F(LK)=F(K)
         F(K)=F0/B(K,K)
         KK=K+1
         F(KK:N)=F(KK:N)-B(KK:N,K)*F(K)
         END DO

      F(N)=F(N)/B(N,N)

         DO K=N1,1,-1
         KK=K+1
            DO J=N,KK,-1
             F(K)=F(K)-B(K,J)*F(J)
            END DO
          END DO

         DO K=N1,1,-1
         MK=MG(K)
         F0=F(MK)
         F(MK)=F(K)
         F(K)=F0
         END DO
100   continue
c
c ERRORS !
      CALL error('GAUSS')
c            
      RETURN
      END
c******************************************************************
c     Метод Ньютона для нахождения стартового решения              *
c     по начальному приближению                                    *
c                                                                  *
      subroutine NEWTON0 (N,NPAR,ndex,limit0,jacob,ifail,exact,
     #                    U,F,U1,U2,DET,STAB,QAPPA,ISTAB)
c
      implicit real*8 (A-H,O-Z)
      implicit integer*4 (i-n)
	include 'number_equations.fi'
	include 'param_model.fi'
c
      REAL*8 U(N+1), F(N+1)
      REAL*8 U1(N+1),  U2(N+1)
      REAL*8 B(NEQ,NEQ)
      REAL*8 V(NEQ), CL(NEQ), Y(NEQ+1)
      COMMON /ADD_R/ S(Nparam)
c---------  Параметры программы  ------------------------------------
      NP1=N+1
      N0=N
      ifail=-1
c
c     Итерации по методу Ньютона
c
      ITER=0
 1001 ITER=ITER+1
c
c   Вычисление элементов матрицы Якоби
c   и правых частей системы. Заполнение вектора
c   производных по параметру
c
        DO I=1,NP1
         if(U(I).lt.U1(I).or.U(I).gt.U2(I)) ifail=14
        END DO
	if(ifail.gt.0) return
c
        S(NPAR)=U(NP1)
        X=S(NPAR)
       CALL PEDERV(N,X,U,ML,MU,B,N0)
       CALL COLUM_Q(N,NPAR,U,CL)
       CALL DIFFUN (N,X,U,F)
      if(ndex.le.N) then 
         V(1:N)=B(1:N,ndex)
         B(1:N,ndex)=CL(1:N)
      else
         V(1:N)=CL(1:N)
	end if
c
c     Вычисление невязок -- C0
c
      CALL GAUSS(N,B,F,IGAUSS,1)
c
      IF(IGAUSS.EQ.1) THEN
        ifail=1 ! Матрица Якоби вырождена !
        RETURN
      END IF
 
       C0=0d0
       C1=0d0
        DO I=1,N
         C0=DMAX1(C0,DABS(F(I)))
         C1=DMAX1(C1,DABS(U(I)))
        END DO
        IF(C1.GT.1d0) C0=C0/C1

      F(NP1)=F(ndex)
      F(ndex)=0d0

      U(1:NP1)=U(1:NP1)-F(1:NP1)
c
c     Проверка на окончание итераций по Ньютону
c
      IF(C0.LT.EXACT) GOTO 200 !!! Выход из Ньютона !

      if(ITER.LT.LIMIT0) GOTO 1001
        ifail=2 !
        RETURN
200   CONTINUE
c     Вычисление производных решения по параметру
c
      F(1:N)=-V(1:N)
c
      CALL GAUSS(N,B,F,IGAUSS,0)
c
      F(NP1)=F(ndex)
      F(ndex)=1d0
c
      if(ISTAB.eq.1) then
       if(STAB.gt.0.d0) then
        X0=S(NPAR)
        S(NPAR)=U(NP1)
        X=S(NPAR)
        Y(1:NP1)=U(1:NP1)	  
        CALL PEDERV(N,X,Y,ML,MU,B,N0)
        CALL STABIL(N,B,STAB,QAPPA)
        CALL CALCDET(N,B,DET)
        S(NPAR)=X0
       end if
      end if
      CALL error('NEWTON0')
c
      RETURN
      END
c--------------------------------------------------------------------------------------
c Calculating the determinant of Jacoby matrix B
c
      SUBROUTINE CALCDET(N,B,DETB)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (i-n)
      DIMENSION B(N,N)
      allocatable P (:)
      allocatable LG(:)
      allocatable MG(:)
c
      allocate ( P(N))
      allocate (LG(N))
      allocate (MG(N))
c-----------------------------------------------------------
      N1=N-1
      DO K=1,N1
       C0=0d0
       DO I=K,N
        DO J=K,N
         IF(DABS(B(I,J)).LT.DABS(C0)) cycle
          C0=B(I,J)
          L0=I
          LG(K)=I
          M0=J
          MG(K)=J
        END DO
       END DO
c
        P(1:N)=B(1:N,M0)
        B(1:N,M0)=B(1:N,K)
        B(1:N,K)=P(1:N)
        P(K:N)=B(L0,K:N)
        B(L0,K:N)=B(K,K:N)
        B(K,K)=C0
        KK=K+1
        DO J=KK,N
         B(K,J)=P(J)/C0
         
	   B(KK:N,J)=B(KK:N,J)-B(KK:N,K)*B(K,J)
        END DO
      END DO
c-----------------------------------------------------------
      BN=B(N,N)
      D0=DSIGN(1.d0,BN)
c      DETB=DLOG(DABS(BN))
      DO I=1,N1
       BI=B(I,I)
       IF(LG(I).NE.I) D0=-D0
       IF(MG(I).NE.I) D0=-D0
       IF(DSIGN(1.d0,BI).LT.0d0) D0=-D0
c       DETB=DETB+DLOG(DABS(BI))
      END DO
c      DETB=D0*DEXP(DETB)
      DETB=D0
      deallocate (P)
      deallocate (LG)
      deallocate (MG)
      RETURN
      END
c*******************************************************************
c     Solving the system of nonlinear equations                    *
c                                                                  *
      subroutine NEWTON(N,NPAR,LIST,ndex,limit,jacob,ipredict,
     *                  hstep,H1,exact,emin,ifail,
     *                  H,U1,U2,U,Y1,Y2,W1,W2,F,DET,STAB,QAPPA) 
c
      implicit real*8 (A-H,O-Z)
      implicit integer*4 (i-n)
	include 'number_equations.fi'
	include 'param_model.fi'
c
      REAL*8 U1(N+1), U2(N+1),U(N+1),H(N+1),F(N+1),Y1(N+1),Y2(N+1), 
     *       W1(N+1), W2(N+1)
      REAL*8 B(NEQ,NEQ) 
      REAL*8 V(NEQ), CL(NEQ), Y(NEQ+1)
      COMMON /ADD_R/   S(Nparam)
c---------  Program constants  ------------------------------------
      NP1=N+1
      N0=N

 1000 ITER=0
      ifail=-1
 1001 ITER=ITER+1
c
c  Проверка на принадлежность решения
c  области определения в текущей позиции
c
      DO I=1,NP1
       IF (U(I).LT.U1(I).OR.U(I).GT.U2(I)) ifail=8
      END DO
      if(ifail.gt.0) goto 14
c
c   Вычисление элементов матрицы Якоби
c   и правых частей системы. Заполнение вектора
c   производных по параметру
c
       S(NPAR)=U(NP1)
       X=S(NPAR) 
       CALL PEDERV(N,X,U,ML,MU,B,N0)
       CALL COLUM_Q(N,NPAR,U,CL)
       CALL DIFFUN (N,X,U,F)
c
      if(ndex.le.N) then
       V(1:N)=B(1:N,ndex)
       B(1:N,ndex)=CL(1:N)	 
      else
       V(1:N)=CL(1:N)
      end if
c
c     Вычисление невязок
c
      CALL GAUSS(N,B,F,IGAUSS,1)
c
      IF(IGAUSS.EQ.1) ifail=6
      if(ifail.gt.0) goto 14

       C0=0d0
       C1=0d0
        DO I=1,N
         C0=DMAX1(C0,DABS(F(I)))
         C1=DMAX1(C1,DABS(U(I)))
        END DO
        IF(C1.GT.1d0) C0=C0/C1

      F(NP1)=F(NDEX)
      F(NDEX)=0d0

      U(1:NP1)=U(1:NP1)-F(1:NP1)
c
c     Проверка на окончание итераций по Ньютону
c
      IF(C0.LT.EXACT) GOTO 200 ! выход из Ньютона  !
      IF(ITER.LT.LIMIT) GOTO 1001
c
c     Деление текущего шага HSTEP пополам
c
   14 HSTEP=5.d-1*HSTEP
      IF(DABS(HSTEP).LT.EMIN*H(NDEX)) THEN
	 if(ifail.lt.0) ifail=7
         RETURN
      END IF
c
      IF(LIST.EQ.1) THEN
        U(1:NP1)=Y2(1:NP1)+HSTEP*W2(1:NP1)
      ELSE
c
c     Вычисление начального приближения в зависимости
c     от порядка аппроксимации 
c
      CALL APPROX (N,ipredict,hstep,H1,U,F,Y1,Y2,W1,W2)
      END IF
      GOTO 1000
200   CONTINUE
c
c     Вычисление производных решения по параметру
c
      F(1:N)=-V(1:N)
      CALL GAUSS(N,B,F,IGAUSS,0)
c
      F(NP1)=F(NDEX)
      F(NDEX)=1d0
c
      if(STAB.gt.0d0) then
       X0=S(NPAR)
       S(NPAR)=U(NP1)
       X=S(NPAR)
       Y(1:NP1)=U(1:NP1)	  
       CALL PEDERV(N,X,Y,ML,MU,B,N0)
       CALL STABIL(N,B,STAB,QAPPA)
       CALL CALCDET(N,B,DET)
       S(NPAR)=X0
      end if
c
      CALL error('NEWTON')
      RETURN
      END
c*********************************************************************  
	Subroutine One_Step(NReturn,Nsys, Nturn,
     *        N,NDEX,IPARAM,Jmatr,NEND,LIMIT,LIMIT0,ifail,
     *        IPREDICT,NPAR,LIST,KBIF,NSOLV,
     *        ISTART, IA, NP1,N0,
     *        HSTEP,ACCEL,EXACT,EMIN,EMAX,H1,
     *        QSOLV,QMIN,QMAX,STAB,QAPPA,EXACT0,
     *        H,
     *        YY, 
     *        U, 
     *        U1,
     *        U2,
     *        Ystep, YQS,
     *        YBIF,
     *        F  ,
     *        Y1 ,
     *        Y2 ,
     *        W1 ,
     *        W2  )

      implicit integer*4 (i-n)
      implicit real*8 (A-H,O-Z)
      include 'number_equations.fi'
      include 'param_model.fi'
c
      REAL*8 YQS (1:Nturn,1:NEQ+2)
      REAL*8 YBIF(1:Nturn,1:NEQ+2)
      REAL*8 Ystep(1:NEQ+4), 
     *       F  (1:NEQ+1),
     *       Y1 (1:NEQ+1),
     *       Y2 (1:NEQ+1),
     *       W1 (1:NEQ+1),
     *       W2 (1:NEQ+1),
     *       H  (1:N+1),
     *       YY (1:N+1), 
     *       U(1:N+1), 
     *       U1(1:N+1),
     *       U2(1:N+1)

	DIMENSION IA(1:Nturn) 
c
      COMMON /ADD_R/ S(Nparam)
c
c------- NReturn = 0 - запись
c------- NReturn = 1 - чтение
c 
	if (NReturn.eq.0) then ! NReturn = 0 - запись
	open(unit=122,file='one_results.bin',FORM='BINARY',
     *     BUFFERED='YES', BLOCKSIZE=204800,
     *     ACCESS='SEQUENTIAL')
	write(122) 
     *        N,NDEX,IPARAM,Jmatr,NEND,LIMIT,LIMIT0,ifail,
     *        IPREDICT,NPAR,LIST,KBIF,NSOLV,
     *        ISTART,IA(1:Nturn),NP1,N0,
     *        HSTEP,ACCEL,EXACT,EMIN,EMAX,H1,
     *        QSOLV,QMIN,QMAX,STAB,QAPPA,EXACT0,
     *        H  (1:NEQ+1),
     *        YY (1:NEQ+1), 
     *        U(1:NEQ+1), 
     *        U1(1:NEQ+1),
     *        U2(1:NEQ+1),
     *        Ystep(1:NEQ+4), YQS(1:Nturn,1:NEQ+2),
     *        YBIF(1:Nturn,1:NEQ+2),
     *        F  (1:NEQ+1),
     *        Y1 (1:NEQ+1),
     *        Y2 (1:NEQ+1),
     *        W1 (1:NEQ+1),
     *        W2 (1:NEQ+1)

 	close(unit=122,STATUS = 'keep')
	else             !NReturn = 1 - чтение
  	open(unit=125,file='one_results.bin', FORM='BINARY',
     *     BUFFERED='YES', BLOCKSIZE=204800,
     *     ACCESS='SEQUENTIAL')
	read(125) 
     *        N,NDEX,IPARAM,Jmatr,NEND,LIMIT,LIMIT0,ifail,
     *        IPREDICT,NPAR,LIST,KBIF,NSOLV,
     *        ISTART,IA(1:Nturn),NP1,N0,
     *        HSTEP,ACCEL,EXACT,EMIN,EMAX,H1,
     *        QSOLV,QMIN,QMAX,STAB,QAPPA,EXACT0,
     *        H  (1:NEQ+1),
     *        YY (1:NEQ+1), 
     *        U(1:NEQ+1), 
     *        U1(1:NEQ+1),
     *        U2(1:NEQ+1),
     *        Ystep(1:NEQ+4), YQS(1:Nturn,1:NEQ+2),
     *        YBIF(1:Nturn,1:NEQ+2),
     *        F  (1:NEQ+1),
     *        Y1 (1:NEQ+1),
     *        Y2 (1:NEQ+1),
     *        W1 (1:NEQ+1),
     *        W2 (1:NEQ+1)

 	close(unit=125,STATUS = 'delete')
	end if
c
	return
	end
c*******************************************************************
c     Solution approximation by linear, квадратичной or cubic      *
c     parabola                                                     *
      subroutine APPROX(N,ipredict,hstep,H1,U,F,Y1,Y2,W1,W2)
c
      implicit real*8 (A-H,O-Z)
      implicit integer*4 (i-n)
c
      REAL*8 U(N+1), F(N+1), Y1(N+1), Y2 (N+1), W1(N+1), W2 (N+1)
c
      NP1=N+1
c
c     Calculating the initial approximation in dependence
c     on the type approximation (IPREDICT=1 (linear),2,3)

      SELECT CASE(ipredict)
      CASE(1)
      U(1:NP1)=Y2(1:NP1)+HSTEP*W2(1:NP1)
      CASE(2) 
      HH=HSTEP/H1
      H0=1.d0+HH
      HY1=HH*(H0-1.d0)
      HY2=H0*(1.d0-HH)
      HW2=H0*HSTEP
      U(1:NP1)=HY1*Y1(1:NP1)+HY2*Y2(1:NP1)+HW2*W2(1:NP1)
      CASE(3)
      HH=HSTEP/H1
      H2=HH*HH
      H3=HH*H2
      HS1=HSTEP*HH
      H0=1.d0+HH
      HY1=3.d0*H2+2.d0*H3
      HY2=1.d0-HY1
      HW1=HS1*H0
      HW2=HSTEP*(1.d0+2.d0*HH+H2)
      U(1:NP1)=HY2*Y2(1:NP1)+HW2*W2(1:NP1)+HY1*Y1(1:NP1)+HW1*W1(1:NP1)
      END SELECT
c
c Errors!
      CALL error('APPROX')
c      
      RETURN
      END
c--------------------------------------------------------------------------------------
      SUBROUTINE BIFUR(N,U,F,Y,W,NDEX,LIST,IBIF,YBIF)
c
      IMPLICIT INTEGER*4 (i-n)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION U(N+1),F(N+1),Y(N+1),W(N+1)
      common /bifur_turn_steps/ Nbif, Nturn, Nsys 
      REAL*8 YBIF(Nturn,N+2)
c
      NP1=N+1
      EPSB=1.D-7
c
      Y0=Y(NP1)
      U0=U(NP1)
      W0=W(NP1)
      F0=F(NP1)

      HH=U(NDEX)-Y(NDEX)

      T1=0.d0
      T2=1.d0
      R1=W0
      R2=F0

   20 IF (T2-T1.LT.EPSB.OR.T2.LT.T1) GOTO 10
      T=0.5d0*(T1+T2)
      T0=1.d0-T
      R0=6.d0*T*T0*(U0-Y0)/HH
      R0=R0+T0*(1.d0-3.d0*T)*W0-T*(2.d0-3.d0*T)*F0
         IF (R0*R2.GT.0d0) GOTO 30
      T1=T
      R1=R0
      GOTO 20
   30 T2=T
      R2=R0
      GOTO 20
   10 CONTINUE
      IBIF=IBIF+1
      TT=T*T
      T00=T0*T0
      T12=(1.d0+2.d0*T)*T00
      T32=(3.d0-2.d0*T)*TT
      HHW=HH*T*T00
      HHF=HH*T0*TT

      YBIF(IBIF,1:NP1)=Y(1:NP1)*T12+U(1:NP1)*T32+
     #                 HHW*W(1:NP1)-HHF*F(1:NP1)
      YBIF(IBIF,NP1+1)=LIST+1
c
c Ошибки
      CALL error('BIFUR')
c 
      RETURN
      END
c--------------------------------------------------------------------------------------
      SUBROUTINE SOLV_Q(N,U,F,Y,W,NDEX,NSOLV,QSOLV,YQS)
      implicit integer*4 (i-n)
      implicit real*8 (A-H,O-Z)
c
      DIMENSION U(N+1),F(N+1),Y(N+1),W(N+1)
      common /bifur_turn_steps/ Nbif, Nturn, Nsys 
      REAL*8 YQS(Nturn,N+2)
c
      NP1=N+1
      EPSQ=1.D-7
c
      Y0=Y(NP1)
      W0=W(NP1)
      U0=U(NP1)
      F0=F(NP1)

      R1=Y0-QSOLV
      R2=U0-QSOLV

      HH=U(NDEX)-Y(NDEX)
c
      T1=0.d0
      T2=1.d0

  20  T=0.5d0*(T1+T2)
      T0=1.d0-T
      T00=T0*T0
      TT=T*T
      T12=(1.d0+2.d0*T)*T00
      T32=(3.d0-2.d0*T)*TT
      HHW=HH*T*T00
      HHF=HH*TT*T0

      R0=Y0*T12+U0*T32
      R0=R0+HHW*W0-HHF*F0-QSOLV

      D1=DSIGN(1d0,R0)
      D2=DSIGN(1d0,R2)
         IF (D1-D2.NE.0d0) THEN
         T1=T
         R1=R0
         ELSE
         T2=T
         R2=R0
         END IF
         IF (T2-T1.LT.EPSQ.OR.T2.LT.T1) GOTO 10
      GOTO 20
   10 NSOLV=NSOLV+1
      YQS(NSOLV,1:N)=Y(1:N)*T12+U(1:N)*T32+HHW*W(1:N)-HHF*F(1:N)
      YQS(NSOLV,NP1)=QSOLV
c
c  ERRORS !
      CALL error('SOLV_Q')
c
      RETURN
      END
c====================================================================
      SUBROUTINE RESTAB(N,LIMSTEP,ifail,STAB,KK,NB,HB,ZB,Ystep)
      implicit INTEGER*4 (i-n)
      implicit real*8 (A-H,O-Z)
      include 'number_equations.fi'
      common /bifur_turn_steps/ Nbif, Nturn, Nsys

      REAL*8 Ystep(N+4)
      REAL*8 Ystep_old(NEQ+4)
      REAL*8 Ystep_new(NEQ+4)
	REAL*8 HB(Nbif*2)
	REAL*8 ZB(Nbif*2)
      character*8 chu
      DIMENSION NB(Nbif*2)
c      
      NP1=N+1
      N0=N
      KK=0
      ifail=-1
c
      OPEN(UNIT=33,FILE='ResTab_old.bin',status='OLD',
     #     access='append', BUFFERED = 'YES',
     #     BLOCKSIZE=20480,FORM='BINARY')
      rewind(33)
	read(33) Ystep_old(:) 
c
         do LIST=2,LIMSTEP
	      if (LIST.eq.LIMSTEP) then
            Ystep_new(:)=Ystep(:)
	      else
	      read(33) Ystep_new(:) 
	      end if

         if(LIST.eq.2) then
          x0=Ystep_old(NP1)
	  IND=NINT(Ystep_old(N+2)) 
          DET0=Ystep_old(N+3)
          QS0=Ystep_old(N+4)

	    if(QS0.lt.0.d0) then
	     chu=' STABLE '
            else
	     chu='UNSTABLE'
            end if
            LL=LIST-1
            write(13,'(5X,I6,e17.8,I4,e12.3,a10)') LL,x0,IND,DET0,chu
         end if

          x1=Ystep_new(N+1)
	  IND =NINT(Ystep_new(N+2)) 
          DET1=Ystep_new(N+3)
          QS1 =Ystep_new(N+4)  

	    if(QS1.lt.0.d0) then
	     chu=' STABLE '
            else
	     chu='UNSTABLE'
            end if
          write(13,'(5X,I6,e17.8,I4,e12.3,a10)') LIST,x1,IND,DET1,chu

          U0=DSIGN(1.d0,QS0)
          D0=DSIGN(1.d0,DET0)
	    
          U1=DSIGN(1.d0,QS1)
          D1=DSIGN(1.d0,DET1)

           if(U0*U1.lt.0.d0) then
            KK=KK+1
	      HB(KK)=x1
	      ZB(KK)=DSIGN(1.d0,D0*D1)
	      NB(KK)=LIST
           if(KK.eq.Nbif*2) then 
            ifail=13
            return
           end if
           end if
c
         Ystep_old(:)=Ystep_new(:)
         QS0=QS1
         DET0=DET1  
         end do
	close(33)
c ERRORS!
      CALL error('RESTAB')
c                    
       RETURN
       END
c--------------------------------------------------------------------
      subroutine reqsolv(NDIM,KSOL,NPAR,STAB,QSOLV,YQS)
       implicit integer*4 (i-n)
       implicit real*8 (A-H,O-Z)
      common /bifur_turn_steps/ Nbif, Nturn, Nsys 
      REAL*8 YQS (Nturn,NDIM+2)
      CHARACTER*8 chu

      NP1=NDIM+1
        write(1,'(/,1x,a,I6)') 'The number of solutions =',KSOL
	write(1,905)'Table 4. Multiplicity solutions and their stability 
     &on st.s.d. '
	write(1,'(5x,a,I4,a,e10.4,/)') 'at the given value of model parame
     &ter p(Npar)=Qsolv, Npar=',NPAR,', Qsolv=',qsolv
        DO J=1,KSOL
          X=YQS(J,NDIM+2) 
          if(X.lt.0d0) then
           write(chu,'(a8)') ' STABLE '
          else
           write(chu,'(a8)') 'UNSTABLE'
          end if
          WRITE(1,'(/,4x,a,I4)') 'Number of solution:',J
          WRITE(1,'(4x,a,a8)') 'Stability: ',chu
          write(1,'(4x,a)') 'Component values:'
          write(1,22) (YQS(J,I),I=1,NDIM)
          write(1,115) ('-',k=1,70)
        END DO
22    FORMAT(2X,100e16.8)
115   FORMAT(70a)
905   FORMAT(/,10x,a,/)
       return
       end 
c
c ----------------------------------------------------------
c
      subroutine sensFM(NDIM,NPAR,FM)
      implicit integer*4 (i-n)
      implicit real*8 (A-H,O-Z)
      include 'number_equations.fi'
      include 'Param_common.fi'
      REAL*8    FM(Nparam,NEQ)	! Модули производных решения по параметрам
      CHARACTER*35    DATLEN,TIMLEN
      COMMON /DATA/   DATLEN,TIMLEN     
C----------------- OUTPUT the table -------------------
      open(unit=1,file='Sensitivity.dat',BUFFERED='YES',
     #     BLOCKSIZE=20480) 
c
      WRITE(1,906) DATLEN
      WRITE(1,906) TIMLEN
	write(1,'(/)')
	write(1,904)'Table 5. Matrix of derivatives of start solution with
     & respect the model parameters (sensitivity matrix)' 
      write(1,'(1x,a,I4,a,I4,a,I4)')
     &                   	'with elements M(i,j)=|dx(j)/dp(i)|, i=1,...
     &,m, m=',Nparam,', j=1,...,n, n=',NDIM,', p(N), N=',NPAR
	write(1,'(/)')
      
	DO J=1,Nparam
	  write(1,999) (FM(J,k),k=1,NDIM)
      END DO  
      
	write(1,'(/,/)')
	write(1,904)'Table 6. Module of maximal derivative in i-th row of 
     &sensitivity matrix,' 
      write(1,'(10x,a,I4)')
     &                   	'max|M(i,j)|=|M(i,j*)|, i=1,...,m, m=',
     &Nparam
	write(1,'(/)')
     
	write(1,115) ('-',k=1,36)     
   
      write(1,'(a)') '   i           j*        |M(i,j*)| '
	write(1,115) ('-',k=1,36)     

	DO J=1,Nparam
	      
	FMAX=0.d0
	imax=1
	DO k=1,NDIM
	if(FMAX.lt.FM(J,k)) then
	FMAX=FM(J,k)
	imax=k
	end if
	END DO
	write(1,'(I4,8X,I4,9X,e10.4)') J,imax,FMAX
	END DO
	write(1,115) ('-',k=1,36)     
      
	write(1,'(/,/)')
	write(1,904)'Table 7. Module of maximal derivative in i-th column  
     &of sensitivity matrix,' 

      write(1,'(10x,a,I4)')
     &                   	'max|M(i,j)|=|M(i*,j)|, j=1,...,n, n=',
     &NDIM
	write(1,'(/)')
     
	write(1,115) ('-',k=1,36)     
   
      write(1,'(a)') '   j           i*        |M(i*,j)| '
	write(1,115) ('-',k=1,36)     

      DO k=1,NDIM
	FMAX=0.d0
	imax=1
	DO J=1,Nparam
	if(FMAX.lt.FM(J,k)) then
	FMAX=FM(J,k)
	imax=J
	end if
	END DO

	write(1,'(I4,8X,I4,9X,e10.4)') k,imax,FMAX

	END DO

	close(1)

115   FORMAT(70a)
904   FORMAT(1x,a,/)
905   FORMAT(/,10x,a,/)
906   FORMAT(A)
999   format(1000(e16.8))
      RETURN
      END
c-----------------------------------------------------------------------------
      subroutine bas_tables (NDIM,KHBP,LIST,NPAR,IHTM,NB,
     #                       STAB,Ystep,HB,ZB)
      implicit integer*4 (i-n)
      implicit real*8 (A-H,O-Z)
      common /bifur_turn_steps/ Nbif, Nturn, Nsys 
      REAL*8   HB(Nbif*2), ZB(Nbif*2), Ystep(NDIM+4)
      DIMENSION NB(Nbif*2)

      CHARACTER*8     chu
      CHARACTER*35    DATLEN,TIMLEN
      COMMON /DATA/   DATLEN,TIMLEN     
C----------------- OUTPUT the table -------------------
      NP1=NDIM+1

!       BASIC TABLES !

	OPEN(UNIT=13,FILE='HTMLresults.dat',
     #     BUFFERED = 'YES',BLOCKSIZE=20480)
      WRITE(13,906) DATLEN
      WRITE(13,906) TIMLEN
      write(13,'(/,a,i7)')' Number of steps, N=',LIST

      WRITE(13,905) '  Table 1. Stationary solutions diagram (st.s.d.))'
      WRITE(13,1112) NPAR 
      if(STAB.gt.0.d0) then 
      WRITE(13,904)
     #	  '    Nstep    Parameter   Index  Sign of DET  Stability'
      CALL RESTAB (NDIM,LIST,ifail,STAB,KHBP,NB,HB,ZB,Ystep)
      if(ifail.eq.13) write(13,906)
     #'    The number of bifurcation points is larger than given one!' 
      else
      WRITE(13,904)
     #	  '    Nstep    Parameter    Index'

      OPEN(UNIT=33,FILE='ResTab_old.bin',status='OLD',
     #     access='append', BUFFERED = 'YES',
     #     BLOCKSIZE=20480,FORM='BINARY')
      rewind(33)
c
         do LL=1,LIST
	    read(33) Ystep(:) 
            IND=NINT(Ystep(NP1+1)) 
            write(13,'(5X,I6,e17.8,I4)') LL,Ystep(NP1),IND
	 end do
         close(unit=33)
      end if 
c
       if(IHTM.lt.1) goto 777 
!                         HTML-файл
	OPEN(UNIT=12,FILE='HTMLresults.html',BUFFERED = 'YES',
     #     BLOCKSIZE=20480)
	WRITE(12,906) '<HTML>'
      WRITE(12,906) '<HEAD>'
      WRITE(12,906) '<TITLE>Stationary solutions diagram</TITLE>'
      WRITE(12,906) '</HEAD>'
      WRITE(12,906) '<BODY BGCOLOR=#E3EFDE>' ! цвет фона документа

      WRITE(12,906) '<BIG>'
      WRITE(12,906) '<P>'

      WRITE(12,906) '<TABLE WIDTH=80%>'
      WRITE(12,906) '<CAPTION ALIGN=TOP>Stationary solutions diagram
     #</CAPTION>'
      WRITE(12,906) '<TR>'

      OPEN(UNIT=33,FILE='ResTab_old.bin',status='OLD',
     #     access='append', BUFFERED = 'YES',
     #     BLOCKSIZE=20480,FORM='BINARY')
      rewind(33)

      rewind(13)
        do i=1,11
        READ(13,*)
	  end do
      IF (STAB.gt.0.d0) THEN

      WRITE(12,906)
     #'<TD>NSTEP</TD><TD>PARAMETER</TD><TD>INDEX</TD><TD>Sign of DET</TD
     #><TD>STABILITY</TD>'

	do K=1,LIST
      read(13,'(5X,I6,e17.8,I4,e12.3,a10)')
     #     l_LIST,z_Ystep_LIST_NP1,
     #     l_IND,z_Ystep_LIST_NP1_3,chu
      WRITE(12,906) '<TR>'                                     

       if(K.lt.10) then
      WRITE(12,'(a,I1,a,I4,a)')'<TD><A HREF="#s0000',K,'">',    
     #l_LIST,'</A></TD>'
        else
       if(K.lt.100) then
      WRITE(12,'(a,I2,a,I4,a)')'<TD><A HREF="#s000',K,'">',    
     #l_LIST,'</A></TD>'
        else
       if(K.lt.1000) then
      WRITE(12,'(a,I3,a,I4,a)')'<TD><A HREF="#s00',K,'">',    
     #l_LIST,'</A></TD>'
        else
       if(K.lt.10000) then
      WRITE(12,'(a,I4,a,I4,a)')'<TD><A HREF="#s0',K,'">',    
     #l_LIST,'</A></TD>'
        else
      WRITE(12,'(a,I5,a,I4,a)')'<TD><A HREF="#s',K,'">',    
     #l_LIST,'</A></TD>'
        end if
        end if
        end if
        end if
      write(12,'(a,e15.8,a,I4,a,e12.3,a10,a)') 
     #      '<TD>',z_Ystep_LIST_NP1,'</TD><TD>',l_IND,
     #      '</TD><TD>',z_Ystep_LIST_NP1_3,'</TD><TD>',chu,'</TD>' 
      WRITE(12,906) '</TR>'   
	end do 
      ELSE
      WRITE(12,906)
     #'<TD>NSTEP</TD><TD>PARAMETER</TD><TD>INDEX</TD>'

	do K=1,LIST
      read(13,'(5X,I6,e17.8,I4)') l_LIST,z_Ystep_LIST_NP1,l_IND
      WRITE(12,906) '<TR>'                                     

       if(K.lt.10) then
      WRITE(12,'(a,I1,a,I4,a)')'<TD><A HREF="#s0000',K,'">',    
     #l_LIST,'</A></TD>'
        else
       if(K.lt.100) then
      WRITE(12,'(a,I2,a,I4,a)')'<TD><A HREF="#s000',K,'">',    
     #l_LIST,'</A></TD>'
        else
       if(K.lt.1000) then
      WRITE(12,'(a,I3,a,I4,a)')'<TD><A HREF="#s00',K,'">',    
     #l_LIST,'</A></TD>'
        else
       if(K.lt.10000) then
      WRITE(12,'(a,I4,a,I4,a)')'<TD><A HREF="#s0',K,'">',    
     #l_LIST,'</A></TD>'
        else
      WRITE(12,'(a,I5,a,I4,a)')'<TD><A HREF="#s',K,'">',    
     #l_LIST,'</A></TD>'
        end if
        end if
        end if
        end if
      write(12,'(a,e15.8,a,I4,a)') 
     #      '<TD>',z_Ystep_LIST_NP1,'</TD><TD>',l_IND,'</TD>'
      WRITE(12,906) '</TR>'   
	end do 
      END IF

      WRITE(12,906) '</TABLE>'                                        
      WRITE(12,906) '<HR SIZE=1 WIDHT=80% ALIGN=LEFT>' ! гориз.линия 

      WRITE(12,906) '</P>'                                           
      WRITE(12,906) '<P>'
      WRITE(12,906) '<TABLE WIDTH=100%>'                             

      rewind(13)
        do i=1,11
        READ(13,*)
	  end do
      IF (STAB.gt.0d0) THEN
      do K=1,LIST
      read(13,'(5X,I6,e17.8,I4,e12.3,a10)')
     #       l_LIST,z_Ystep_LIST_NP1,
     #       l_IND,z_Ystep_LIST_NP1_3,chu
      WRITE(12,906) '<TR>'                                             

       if(K.lt.10) then
      write(12,'(a,I1,a,I4,a,e17.8,a)') 
     #  '<TD BGCOLOR=#E3BC6C><A NAME=#s0000', 
     #	K,'>',l_LIST,'</A></TD><TD>',z_Ystep_LIST_NP1,'</TD>'
       else 
       if(K.lt.100) then
      write(12,'(a,I2,a,I4,a,e17.8,a)') 
     #  '<TD BGCOLOR=#E3BC6C><A NAME=#s000', 
     #	K,'>',l_LIST,'</A></TD><TD>',z_Ystep_LIST_NP1,'</TD>'
       else 
       if(K.lt.1000) then
      write(12,'(a,I3,a,I4,a,e17.8,a)') 
     #  '<TD BGCOLOR=#E3BC6C><A NAME=#s00', 
     #	K,'>',l_LIST,'</A></TD><TD>',z_Ystep_LIST_NP1,'</TD>'
       else 
       if(K.lt.10000) then
      write(12,'(a,I4,a,I4,a,e17.8,a)') 
     #  '<TD BGCOLOR=#E3BC6C><A NAME=#s0', 
     #	K,'>',l_LIST,'</A></TD><TD>',z_Ystep_LIST_NP1,'</TD>'
       else 
      write(12,'(a,I5,a,I4,a,e17.8,a)') 
     #  '<TD BGCOLOR=#E3BC6C><A NAME=#s', 
     #	K,'>',l_LIST,'</A></TD><TD>',z_Ystep_LIST_NP1,'</TD>'
        end if
        end if
        end if
        end if
      read(33) Ystep(:) 
	DO ii=1,NDIM
	write(12,'(a,e15.8,a)') '<TD>',Ystep(ii),'</TD>'            
	END DO
	WRITE(12,906) '</TR>'
	end do
        ELSE
      do K=1,LIST
      read(13,'(5X,I6,e17.8,I4)')
     #       l_LIST,z_Ystep_LIST_NP1,l_IND
      WRITE(12,906) '<TR>'                                             

       if(K.lt.10) then
      write(12,'(a,I1,a,I4,a,e17.8,a)') 
     #  '<TD BGCOLOR=#E3BC6C><A NAME=#s0000', 
     #	K,'>',l_LIST,'</A></TD><TD>',z_Ystep_LIST_NP1,'</TD>'
       else 
       if(K.lt.100) then
      write(12,'(a,I2,a,I4,a,e17.8,a)') 
     #  '<TD BGCOLOR=#E3BC6C><A NAME=#s000', 
     #	K,'>',l_LIST,'</A></TD><TD>',z_Ystep_LIST_NP1,'</TD>'
       else 
       if(K.lt.1000) then
      write(12,'(a,I3,a,I4,a,e17.8,a)') 
     #  '<TD BGCOLOR=#E3BC6C><A NAME=#s00', 
     #	K,'>',l_LIST,'</A></TD><TD>',z_Ystep_LIST_NP1,'</TD>'
       else 
       if(K.lt.10000) then
      write(12,'(a,I4,a,I4,a,e17.8,a)') 
     #  '<TD BGCOLOR=#E3BC6C><A NAME=#s0', 
     #	K,'>',l_LIST,'</A></TD><TD>',z_Ystep_LIST_NP1,'</TD>'
       else 
      write(12,'(a,I5,a,I4,a,e17.8,a)') 
     #  '<TD BGCOLOR=#E3BC6C><A NAME=#s', 
     #	K,'>',l_LIST,'</A></TD><TD>',z_Ystep_LIST_NP1,'</TD>'
        end if
        end if
        end if
        end if
      read(33) Ystep(:) 
	DO ii=1,NDIM
	write(12,'(a,e15.8,a)') '<TD>',Ystep(ii),'</TD>'            
	END DO
	WRITE(12,906) '</TR>'
	end do
        END IF
c---------------------------------------
        close(unit=33)

	DO JJ=1,500                                                        
         WRITE (12,906) '<TR><TD> </TD></TR>'
	END DO

      WRITE(12,906) '</TABLE>'
      WRITE(12,906) '</P>'

      WRITE(12,906) '</BIG>'
	WRITE(12,906) '</BODY>'
      WRITE(12,906) '</HTML>'

	CLOSE(UNIT=12)
777     CONTINUE
        CLOSE(UNIT=13)

904   FORMAT(1x,a,/)
905   FORMAT(/,10x,a,/)
906   FORMAT(A)
1112  FORMAT(5X,'Results of the model parameter continuati
     #on p(N), N=',I4,/)
      RETURN
      END 
c!!!! END PROGRAM  RISOLVE, version from 23.03.2007
