Program to numerically solve 3 ODE's

Help Homework
26

It looks good, but I have inserted some suggestions and questions as comments prefaced by !!.

Program AE6009
	implicit none
!
	integer, parameter :: dp = kind(1d0)
	integer, parameter :: n = 3						!Order of ODE's
	real(dp) :: a, b, c, g							!Variables for shooting method
	real(dp) :: eta, deta, etastart, etastop		!Variables for Runge-cutta
	real(dp), dimension(n) :: y, k1, k2, k3, k4
!	
	Write (6,*)'' !! use * instead of 6 or output_unit from iso_fortran_env
	write (6,*) 'Enter a guess for f"(0) < 0'
	Read (*,*) a
	Write (6,*) 'Enter a guess for f"(0) > 0'
	Read (*,*) b
!	
	g = 1.0d0
	do while (g > 1d-10 .or. g < -1d-10)  
!
	C=0.5*(a+b) !! why is the do while block not indented? Also, use 0.5_dp instead of 0.5
	etastart = 0.0d0 !! this group of statements is repeated below. It could be made 
                     !! an internal subroutine
	etastop = 4.8d0		! eta at infinity
	deta = 0.01d0
	y(1) = 0.0d0		! Initial Condition f(0) = 0
	y(2) = 0.0d0		! Initial Condition f'(0) = 0
	y(3) = c			! Initial Condition guess for f''(0)
	eta = etastart
!
		do while(eta < etastop)
			call Runge ()
			if (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-10) then
			end if !! should there be a statement before this?
		end do
!
	g = y(2)-1d0
!
		if (g > 0d0)then !! put a space before then for readability
			a = c
			else !! should have same indentation as if
			b = c
		end if
	end do
	write (6,*)''
	Write (6,*) 'The value for f"(0) that satisfies fprime(inf) = 1'
	Write (6,'(a,f10.5)') 'f"(0) =' ,c
	Write (6,*)''
!
! 	Now calculate the solutions to the ODE's
!
	etastart = 0.0d0
	etastop = 4.8d0		! eta at infinity
	deta = 0.01d0
	y(1) = 0.0d0		! Initial Condition f(0) = 0
	y(2) = 0.0d0		! Initial Condition f'(0) = 0
	y(3) = c			! Initial Condition f"(0) = c
	eta = etastart
	Write (6, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
!
	do while(eta < etastop)
			call Runge () !! should have same indentation as line below
		if (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-12) then
			Write (6, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
		end if
	end do
!
!	Program subroutines and sub functions 
!
Contains
	Subroutine Runge ()
		k1 = deta*f(eta,y)						!Runge-cutta 4th order formulas
		k2 = deta*f(eta+deta/2, y+k1/2)
		k3 = deta*f(eta+deta/2, y+k2/2)
		k4 = deta*f(eta+deta, y+k3)
		y = y+((1d0/6)*(k1+(2*k2)+(2*k3)+k4))
		eta = eta+deta
	end subroutine
!
	Function f(eta,y)
		real(dp), intent(in) :: eta, y(n)
		real (dp) :: f(n)
		f = [y(2), y(3), y(2)*y(2)-1d0-y(1)*y(3)] ! Hiemenz Equation
	end Function
!
end Program AE6009
1 Like