!
! Program to solve 2f''' + f f'' = 0
! Using RK-4. For midterm exam, MAE 185, spring 2006
!
! Initial conditions at t=0 are f=0, f'=0, f''=0.322
!
! Stop iteration when f'' value becomes less than machine epsilon. 
!
! by Nasser Abbasi 

PROGRAM mid_term_spring_2006_MAE185
IMPLICIT NONE

!*************************
!*  DATA DECLARATION 
!**************************
REAL *8 :: machine_eps = EPSILON(1.0d0),x1=0.0,x2=0.0,x3=0.322,t=0.0,h,h2
REAL *8 :: x1_prime,x2_prime,x3_prime
INTEGER :: step=0,PRINT_PER_STEP

PARAMETER(h=0.000001, h2=h/2., PRINT_PER_STEP=10000)

!*************************
!*  EXECUTION PART
!**************************

   DO WHILE( x3 > machine_eps )
       
      CALL make_one_step
      t = t + h

      step = step + 1
      IF ( MOD(step,PRINT_PER_STEP) ==0 ) THEN   
          WRITE( UNIT=*,FMT='(4(F24.18,1X))' ) t, x1, x2, x3
          step = 0
      END IF

   END DO


CONTAINS

   SUBROUTINE make_one_step
      ! use RK-4 with interleaving.

      IMPLICIT NONE

      REAL *8 :: k1_x1, k1_x2, k1_x3
      REAL *8 :: k2_x1, k2_x2, k2_x3
      REAL *8 :: k3_x1, k3_x2, k3_x3
      REAL *8 :: k4_x1, k4_x2, k4_x3

      k1_x1 = x1_prime( t,x1,x2,x3 )
      k1_x2 = x2_prime( t,x1,x2,x3 )
      k1_x3 = x3_prime( t,x1,x2,x3 )

      k2_x1 = x1_prime( t+h2 , x1+ h2*k1_x1, x2+ h2*k1_x2, x3+ h2*k1_x3 )
      k2_x2 = x2_prime( t+h2 , x1+ h2*k1_x1, x2+ h2*k1_x2, x3+ h2*k1_x3 )
      k2_x3 = x3_prime( t+h2 , x1+ h2*k1_x1, x2+ h2*k1_x2, x3+ h2*k1_x3 )

      k3_x1 = x1_prime( t+h2 , x1+ h2*k2_x1, x2+ h2*k2_x2, x3+ h2*k2_x3 )
      k3_x2 = x2_prime( t+h2 , x1+ h2*k2_x1, x2+ h2*k2_x2, x3+ h2*k2_x3 )
      k3_x3 = x3_prime( t+h2 , x1+ h2*k2_x1, x2+ h2*k2_x2, x3+ h2*k2_x3 )

      k4_x1 = x1_prime( t+h  , x1+ h*k3_x1,  x2+ h*k3_x2,  x3+ h*k3_x3 )
      k4_x2 = x2_prime( t+h  , x1+ h*k3_x1,  x2+ h*k3_x2,  x3+ h*k3_x3 )
      k4_x3 = x3_prime( t+h  , x1+ h*k3_x1,  x2+ h*k3_x2,  x3+ h*k3_x3 )

      x1 = x1 + h/6.0 * (k1_x1 + 2*k2_x1 + 2*k3_x1 + k4_x1)
      x2 = x2 + h/6.0 * (k1_x2 + 2*k2_x2 + 2*k3_x2 + k4_x2)

      x3 = x3 + h/6.0 * (k1_x3 + 2*k2_x3 + 2*k3_x3 + k4_x3)

   END SUBROUTINE

end PROGRAM   

!***************
!*
!*
!***************
real *8 function x1_prime(t,x1,x2,x3)
IMPLICIT NONE

REAL *8 :: t,x1,x2,x3

   x1_prime=x2

end function
!***************
!*
!*
!***************
real *8 function x2_prime(t,x1,x2,x3)
IMPLICIT NONE

REAL *8 :: t,x1,x2,x3

   x2_prime=x3

end function
!***************
!*
!*
!***************
real *8 function x3_prime(t,x1,x2,x3)
IMPLICIT NONE

REAL *8 :: t,x1,x2,x3

  x3_prime= -(dble(1/2.))*x1*x3
  
end function