!      program rungekinc

	external fcn
      open(1,file='inpaula.txt')
	open(2,file='outnum.txt')
	open(3,file='outteo.txt')
	open(4,file='outexemplo.txt')
!
!	integrador de EDO's com passo fixo - RK 4a ordem
!
	read(1,*)n,tau0
	write(*,*)'n=',n,'tau0=',tau0
	read(1,*) tend,dtau,teta0
	write(*,*)'tend=',tend,'dtau=',dtau,'teta0=',teta0
!
!	initial values
!
	erro=0.
	teta=teta0
	tetateo=teta0
	time=tau0
	k=0 
	write(*,*) '     Tabela de resultados'
	write(*,*)'---------------------------------------------------------------------'
	write(*,*)'  Passo Nr     t          x(t)          xteo(t)          Erro (ETL)'
	write(*,*)'---------------------------------------------------------------------'
	write(*,*)k,time,teta,tetateo,erro

	write(4,*) '     Tabela de resultados'
	write(4,*)'---------------------------------------------------------------------'
	write(4,*)'  Passo Nr     t          x(t)          xteo(t)          Erro (ETL)'
	write(4,*)'---------------------------------------------------------------------'
	write(4,*)k,time,teta,tetateo,erro
	write(2,*)time,teta
	write(3,*)time,teta
 50	k=k+1 
	th=time+0.5*dtau 
	f1=dtau*fcn(time,teta)
	f2=dtau*fcn(th,teta+0.5*f1)
	f3=dtau*fcn(th,teta+0.5*f2)
	f4=dtau*fcn(time+dtau,teta+f3)
	teta=teta+(f1+2.*f2+2.*f3+f4)/6.
	time=time+dtau
     	if (time.le.tend) then
			write(*,*)k,time,teta
			write(2,*)time,teta
		!	tetateo=exp(-time**2/2)  
		!	tetateo=time/(0.5+log(time))
			erro=abs(teta-tetateo)
		!	write(3,*)time,tetateo
			write(4,*)k,time,teta,tetateo,erro
			goto 50 
		else
			goto 300
		endif

300 continue
	close(2)
	close(3)
	close(4)


  call system('notepad outnum.txt') 

  call system('notepad outteo.txt') 

  call system('notepad outexemplo.txt') 

  call system('wgnuplot dados.gnu') ! gráfico

	stop
	end
!**************************************
	function fcn(t,f)
	! aux=-t*f  
!	aux=(t*f-f**2)/t/t
!------3TE---- 2017-2 ----- Noite -----------
!
	rho=1.2
	vol=0.01
	cv=714.
	rk=0.026
	rl1=0.1
	rl2=0.2
	aa=1.
	ta=1000.

	aux1=1./(rho*vol*cv)
	aux2=-rk/((rl1+rl2)/2)

	aux=aux1*(aux2*aa*(f-ta))

	fcn=aux
	return
	end
!***************************************





