!      program rungekinc

      implicit real *8 (a-h,o-z)

	parameter (nelmax=100)
	dimension tp(nelmax)
	external fcn,rkqc
      common /const/ ht2 
  	  common /adap/ alfa,p,tol
  	  common /passo/ h

      open(1,file='inpaula.txt')
	open(2,file='x1.txt')
	open(4,file='x2.txt')
	open(8,file='h.txt')
	open(9,file='x-all.txt')
!
!	integrador de EDO's com passo adaptativo - RK 4a/5a ordem
!	integrador de EDO's com passo fixo - RK 4a ordem  e For. Euler
!
	read(1,*)iflag
	write(*,*)'iflag=',iflag
	read(1,*)n,tau0
	write(*,*)'n=',n,'tau0=',tau0
	read(1,*) tend,dtau,teta0
	write(*,*)'tend=',tend,'dtau=',dtau,'teta0=',teta0
	read(1,*) tol
	write(*,*)'tol=',tol
	read(1,*) alfa
	write(*,*)'alfa=',alfa
	read(1,*) p
	write(*,*)'p=',p
	read(1,*) ht2
	write(*,*)'ht2=',ht2
	read(1,*) dtab
	write(*,*)'dtab=',dtab

!
!	initial values
!
	time=tau0
	do i=1,n
		tp(i)=teta0
	enddo
	teta=teta0
	k=0 
	write(9,*) '     Table of results'
	write(9,*)'-----------------------------'
	write(9,*)'  Passo Nr     time       T(i)    '
	write(9,*)'-----------------------------'
	write(2,*)time,tp(1)
	write(4,*)time,tp(2)

!	ht2=dtau/10
	h=ht2

!   flag para produzir tabela numérica
	aux=dtab

	write(8,*)time,h
	write(9,1001)k,time,(tp(l),l=1,n),h

1001 format(i5,1x,f12.6,2f12.6,f12.6)
!
!	beginning of time loop
!
	iflagfim=0

 50	k=k+1 
	tendi=min(time+dtau,tend)
      	write(*,*)'-------------time=',tendi
	if(iflag.eq.0) then	! RK-adaptativo
     call odeint(tp,n,time,tendi,1.e-6,ht2,1.e-20,id1,id2,nelmax,fcn,rkqc)
	endif
	if(iflag.eq.1) then	 ! RK-passo fixo
	      kk=0 
 500	  kk=kk+1 
		time=min(time+ht2,tendi) 
		call  rk4ord(tp,n,time,ht2,fcn,nelmax)	
		if (time.lt.tendi) goto 500 
	endif
	if(iflag.eq.2) then	  ! Forward Euler
	    call fore(n,fcn,time,tp,tendi,nelmax)
	endif
	if(iflag.eq.3) then	  ! RK-adaptativo-alternativo
	    call rk4feadap(n,fcn,time,tp,tendi,nelmax)
		write(8,*)tendi,h
	endif

	write(*,*)k,tendi,(tp(l),l=1,n)
	write(2,*)tendi,tp(1)
	write(4,*)tendi,tp(2)

	if(tendi.gt.aux) then
		write(9,1001)k,tendi,(tp(l),l=1,n),h
		aux=aux+dtab
	endif

  !     	if (tendi.lt.tend) then
     	if (tp(1).lt.1.d0) then
			time=tendi
			goto 50 
		else
			write(9,1001)k,tendi,(tp(l),l=1,n),h
		endif
	close(2)
	close(4)
	close(8)
	close(9)

call system('notepad x1.txt')       ! listagem 1a variavel
call system('notepad x2.txt')       ! listagem 2a variavel
 call system('notepad h.txt')       ! listagem do passo h-variavel no metodo alternativo

 call system('notepad x-all.txt')       ! listagem de tudo

call system('wgnuplot dados.gnu')  ! gráfico	geral
	stop
	end
!**************************************
	subroutine fcn(n,t,fi,f,nelmax)
      implicit real *8 (a-h,o-z)
	dimension fi(nelmax),f(nelmax)
	!	f(1)=(t*fi(1)-fi(1)**2)/t/t+fi(2)
	!	f(2)=fi(2)-fi(1)
	! Exemplo 7.12 do livro - queda de bloco de gelo
		g=1.d4
		b=4
	 	f(1)=fi(2)
		rlam=2*t/(g+sqrt(g**2+4*t))   
		f(2)=1.-(1.-fi(1))*b*fi(2)/rlam
	return
	end
!***************************************
!234567890123456789012345678901234567890123456789012345678901234567890 
      subroutine fore(n,fcn,time,fi,tend,nelmax) 
! 
      implicit real *8 (a-h,o-z)
      parameter (nd1=100) 
      dimension fi(nelmax),f(nd1)
      common /const/ ht2 
      external fcn 
      k=0 
 50   k=k+1 
      time=min(time+ht2,tend) 
      call fcn(n,time,fi,f,nelmax) 
      do 100 i=1,n 
         fi(i)=fi(i)+ht2*f(i) 
 100  continue 
      if (time.lt.tend) goto 50 

  !	write(8,*)'h = ',h

      return 
      end 
!------------------------------------------------------------------ 





