!      program blasius
use msflib !  systemqq

logical calling
!	parameter (nelmax=100)
!	dimension tp(nelmax)
!	external fcn,rkqc
	common /param/ eta0,deta,etaend,f10,f20
	common /noeq/ n					

	external ff,tf

 !     common /const/ ht2 
      open(1,file='inpclass.txt')
	open(2,file='f1.txt')
	open(3,file='f2.txt')
	open(4,file='f3.txt')
!	open(5,file='f2.txt')
!	open(6,file='f3.txt')

!	Solution of a Boundary Value Problem - BVP
!
!	data input
!
	read(1,*)n
	write(*,*)'n=',n
	read(1,*)eta0
	write(*,*)'eta0=',eta0
	read(1,*)etaend
	write(*,*)'etaend=',etaend
	read(1,*)deta
	write(*,*)'deta=',deta
	read(1,*)f10
	write(*,*)'f10=',f10
	read(1,*)f20
	write(*,*)'f20=',f20
	read(1,*)f30a
	write(*,*)'f30a=',f30a
	read(1,*)f30b
	write(*,*)'f30b=',f30b
	read(1,*)maxit
	write(*,*)'maxit=',maxit
	read(1,*)tol
	write(*,*)'tol=',tol

	! Blasius solution

	call secante(f30a,f30b,maxit,tol,c,ff)

	! Energy equation

! read f1.txt storing the f values in a vector f1(eta)

!	call secante(t20a,t20b,maxit,tol,c,tf)

	write(*,*) 'f30=',c

	close(2)
	close(3)
	close(4)

calling = systemqq('notepad f1.txt')       ! list of data

calling = systemqq('notepad f2.txt')       ! list of data
calling = systemqq('notepad f3.txt')       ! list of data

chamada = systemqq('wgnuplot data.gnu')  ! gráfico
	stop
	end
!**************************************
	subroutine fcn(n,t,fi,f,nelmax)
	dimension fi(nelmax),f(nelmax)
!
	f(1)=fi(2)
	f(2)=fi(3)
	f(3)=-fi(1)*fi(3)/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 
!------------------------------------------------------------------ 
	function ff(x)
	parameter (nmax=100)
	dimension tp(nmax)
	external fcn,rkqc
	common /noeq/ n
	common /flag/ iflag

	common /param/ eta0,deta,etaend,f10,f20
!
!	initial values
!
	eta=eta0
	tp(1)=f10
	tp(2)=f20
	tp(3)=x
!
	k=0 
	if (iflag.eq.1) then

	write(2,*)eta,tp(1)
	write(3,*)eta,tp(2)
	write(4,*)eta,tp(3)
	endif
!
	ht2=deta/1
	nd=nmax
!
!	beginning of time loop
!
 50	k=k+1 
	tendi=eta+deta
!      	write(*,*)'-------------eta=',tendi
!
!     call odeint(tp,n,time,tendi,1.e-6,ht2,1.e-20,id1,id2,nelmax,fcn,rkqc)
	call  rk4ord(tp,n,eta,deta,fcn,nd)
!    call fore(n,fcn,time,tp,tendi,nmax)
	if (iflag.eq.1) then
!	write(*,*)k,tendi,(tp(l),l=1,n)
	write(2,*)tendi,tp(1)
	write(3,*)tendi,tp(2)
	write(4,*)tendi,tp(3)
	endif
!
     	if (tendi.lt.etaend) then
			eta=tendi
			goto 50 
		endif

	aux=tp(2)-1.
	ff=aux


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


