*********************************************************************
*****  Calculation of Velocity Boundary Layer over a flat plate *****              
**********************************************************************
     
	
	real A,kf,kY,kZ
	dimension eta(500),Y(500),Z(500),kf(500),kY(500),kZ(500),f(500),
     t	q(500),err(500)
      
	open (unit=1,file='velocity_bl.dat')
	     

*********************************************************************   
***** Equations for Velocity boundary layers over a flat plate  *****
*****                 f'=Y                                      ***** 
*****                 Y'=Z                                      *****
*****                 Z'=-(1./2.)fZ                             *****
*****  BC's: @eta=0,f=f'=0,Z=A(guessed) ; @eta=INF.,Y=1         *****
***** A is the guessed value for Z at eta=0                     *****
***** Runge kuuta method is used to fing f,Y,Z                  *****
*********************************************************************
 
	 A=1.

	 i=0
       h=0.05
   1   Z(1)=A
       i=i+1
   	      
		  
           Y(1)=0 
	     f(1)=0
	     eta(1)=0

      do 2  n=1,150
    	   kf(1)=Y(n)
   	   kY(1)=Z(n)
	   kZ(1)=-(1./2.)*Z(n)*f(n)
	        
	   kf(2)=(Y(n)+(h/2.)*kY(1))
	   kY(2)=(Z(n)+(h/2.)*kZ(1))
	   kZ(2)=-(1./2.)*(Z(n)+(h/2.)*kZ(1))*(f(n)+(h/2.)*kf(1))
	   
	   kf(3)=(Y(n)+(h/2.)*kY(2))
	   kY(3)=(Z(n)+(h/2.)*kZ(2))
	   kZ(3)=-(1./2.)*(Z(n)+(h/2.)*kZ(2))*(f(n)+(h/2.)*kf(2))
	   
	   kf(4)=(Y(n)+h*kY(3))
	   kY(4)=(Z(n)+h*kZ(3))
	   kZ(4)=-(1./2.)*(Z(n)+h*kZ(3))*(f(n)+h*kf(3))
	   
	   eta(n+1)=eta(n)+h
	   
	   f(n+1)=f(n)+h*(kf(1)+2*kf(2)+2*kf(3)+kf(4))/6
	   Y(n+1)=Y(n)+h*(kY(1)+2*kY(2)+2*kY(3)+kY(4))/6
	   Z(n+1)=Z(n)+h*(kZ(1)+2*kZ(2)+2*kZ(3)+kZ(4))/6
        
    2  continue

	

	     q(i)=A
	     err(i)=(1.-Y(91))

***************************************************************************
*****     We Are using Newton Rapson mathod to find A_impproved       *****
*****              A_improved=A-(Y(A)-1)/(dY/dA)                      ***** 
*****     Here,      Y(A)-1=err(i)                                    *****
*****                dY/dA=err(i)-err(i-1)/q(i)-q(i-1)                *****
***************************************************************************
	
	if (i.eq.1) then
	      A=A+0.05
	      goto 1
	else

	     if (i.eq.3) then
	            if (err(2).lt.err(1)) then
                  A_improved=q(2)-((q(3)-q(2))*err(2)/(err(3)-err(2)))
	            else                 
	            A_improved=q(1)-((q(3)-q(1))*err(1)/(err(3)-err(1)))
	            endif

	     else
	     A_improved=q(i-1)-((q(i)-q(i-1))*err(i-1)/(err(i)-err(i-1)))
	     endif
	
	     
	endif
	
	
	     if (abs(err(i)).gt.0.0001) then 
	     A =A_improved
           goto 1
	     else

	 write (1,10)
  10	 format(1x,'SOLUTION OF HYDRODYNAMIC BOUNDARY LAYER FOR FLOW OVER 
     d FLAT PlATE'/)
   
      

	write(1,20)
	write (1,30)
  30  format (4x,'ETA',6x,'[F]',6x,'[FPRIME]',4x,'[F"]')
         
	write (1,20)
  20  format ('____________________________________________________')
	write(1,40)(eta(x),f(x),Y(x),Z(x),x=1,150)	
     
  40  format(f10.6,f10.6,1x,f10.6,2x,f8.6)
	
	     endif
      stop
	end