!      program amortecimento
! use msflib ! biblioteca que contem o comando systemqq

! logical chamada

	implicit real *8 (a-h,o-z)
	parameter (nelmax=100)
	dimension tp(nelmax)
	external fcn,rkqc
      common /const/ ht2 
      open(1,file='inpaula.txt')
	open(2,file='x1.txt')
	open(3,file='x-all.txt')
	open(4,file='x2.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
	write(*,*)'tend=',tend,'dtau=',dtau
!
!	initial values
!
		read(1,*)(tp(l),l=1,n)
		write(*,*)(tp(l),l=1,n)

!	initial time

	time=tau0

	k=0 
	write(3,*) '     Table of results'
	write(3,*)'-----------------------------'
	write(3,*)'  Passo Nr     time       T(i)    '
	write(3,*)'-----------------------------'
	write(3,*)k,time,(tp(l),l=1,n)
	write(2,*)time,tp(1)
	write(4,*)time,tp(2)
	ht2=dtau/100
!
!	beginning of time loop
!
 50	k=k+1 
	tendi=time+dtau
      	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
	write(3,*)k,tendi,(tp(l),l=1,n)
	write(2,*)tendi,tp(1)
	write(4,*)tendi,tp(2)
     	if (tendi.lt.tend) then
			time=tendi
			goto 50 
		endif
	close(2)
	close(3)
	close(4)

call system('notepad x-all.txt')       ! listagem de todas as variaveis
! chamada = systemqq('notepad x2.txt')       ! listagem 2a variavel


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.8 do livro - sistema sub-, critico e superamortecido
		zeta=1.4
		w_n=0.1
		f(1)=fi(2)
		f(2)=w_n**2*(-fi(1)-2*zeta*fi(2)/w_n)
	! Exemplo 7.9 do livro - presa-predador
	!	a=1.
	!	b=2.
	!	c=1.
	!	d=1.
	!	f(1)=-c*fi(1)+d*fi(1)*fi(2)
	!	f(2)=a*fi(2)-b*fi(1)*fi(2)

	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 
      return 
      end 
!------------------------------------------------------------------ 





