Program shooting
!---------------------------------------------------------------------
! Solver for 1D boundary-value problem second-order ODE 
! Method: unilizes the shooting method based on the method of secants
! (calls 4th-order Runge-Kutta to solve the initial value problem)
! Demo version for d2y/dx2 = f(x,y,dy/dx) equation
! AG: Last revision October 2009
!----------------------------------------------------------------------

implicit none

integer, parameter :: n=21                     ! number of base points
double precision, external :: f
double precision, dimension(1:n) :: x, y, dy   ! x, y, y'
integer i

!open (unit=7,file="shoot2.dat")

! boundary values
x(1) =  0.0
y(1) =  0.0
x(n) =  1.0
y(n) =  1.0

! assumptions for y'(1) - use dy(1) and dy(2) here only as a storage
dy(1) =  2.0
dy(2) =  8.0

call shoot2(f,x,y,dy,n)

write(*,100)
do i=1,n
  write (*,101) x(i), y(i), dy(i)
end do

100 format(5x,'x',11x,'y',11x,'dy')
101	format(3(1pe12.4))
stop
end program shooting

  subroutine shoot2(f,x,y,dy,n)
!======================================================================
! Solution of the boundary-value second-order 1D ODE
! d2y/dx2 = f(x,y,dy/dx) with Dirichlet boundary conditions
! y(xmin) = ..., and y(xmax) = ...
! Method: unilizes the shooting method based on the method of secants
! (calls 4th-order Runge-Kutta to solve the initial value problem)
! written by:  Alex Godunov (last revision - 9 October 2009)
!----------------------------------------------------------------------
! input ...
!  f(x,y,dy)  - function d2y/dx2 (supplied by a user)
!  x(1), x(n) - boundary points
!  y(1), y(n) - boundary values (Dirichlet boundary conditions)
! dy(1),dy(2) - two guesses for y'(x(1))
!
! output ...
!  y(i) and dy(i) solutions at points x(i) (i=1,...,n)
! note: dy corresponds to y' (the first derivative)
!======================================================================
implicit none
integer, parameter :: it=21                    ! max numer of iterations
double precision, parameter :: eps = 0.000001  ! target tolerance
integer n
double precision, external:: f
double precision, dimension(1:n) :: x, y, dy
integer i,j
double precision dx, yn
double precision, dimension(1:it):: g, c

! g(it) - array of dy(it) values [iterative improvement]
! c(it) - array of y(n) values [iterative solutions for g(it)]
! first guesses for g(it)
g(1) = dy(1)
g(2) = dy(2)

! remember the second boundary condition 
! since y(n) recalculated for each new g(i)
yn = y(n)

! generate base points x(i) from x(1), x(n) and n
dx = (x(n)-x(1))/float(n-1)
do i=2,n
  x(i) = x(i-1)+dx
end do

! shooting iterations (for the first two - we use assumed values of dy(1))
do j=1,it
  dy(1) = g(j)
  call rk4_2d(f,x,y,dy,n)      ! solves initial value ODE on n-base points
  c(j) = y(n)
  if(abs(yn-c(j)) <= eps) exit
  if(j.ge.2) g(j+1)=g(j)-(c(j)-yn)*(g(j)-g(j-1))/(c(j)-c(j-1)) ! secant method
end do
end subroutine shoot2

  function f(x,y,dy)
!------------------------------------------
! the second derivative - use original ODE
! d2y/dx2 = f(x,y,dy)
!------------------------------------------
  implicit none
  double precision f, x, y, dy
    f = 2.0*y*y
  end function f

  subroutine rk4_2d(f,x,y,dy,n)
!================================================================
! Solution of the second-order 1D ODE (Initial-value problem)
! method:  Runge-Kutta 4th-order
! written by: Alex Godunov (last revision October 2009)
!----------------------------------------------------------------
! input ...
!  f(x,y,dy)- function from d2y/dx2=f(x,y,dy) (supplied by a user)
!  x(i)  - array of x base point (n points)
!  y(1)  - initial value for y(1)
! dy(1) - initial value for y'(1)
!
! output ...
!  y(i)  - solutions for y  in n points
! dy(i)  - solutions for y' in n points 
!=================================================
implicit none
double precision f
integer n
double precision, dimension(1:n) :: x, y, dy
integer i
double precision h,k11,k12,k21,k22,k31,k32,k41,k42

do i=2,n
   h   = x(i)-x(i-1)
   k11 = h*dy(i-1)
   k12 = h*f(x(i-1),y(i-1),dy(i-1))
   k21 = h*(dy(i-1)+k12/2.0)
   k22 = h*f(x(i-1)+h/2.0, y(i-1)+k11/2.0, dy(i-1)+k12/2.0)
   k31 = h*(dy(i-1)+k22/2.0)
   k32 = h*f(x(i-1)+h/2.0, y(i-1)+k21/2.0, dy(i-1)+k22/2.0)
   k41 = h*(dy(i-1)+k32)
   k42 = h*f(x(i-1)+h,y(i-1)+k31,dy(i-1)+k32)

   y(i) = y(i-1) + (k11 + 2.0*(k21+k31) + k41)/6.0
  dy(i) = dy(i-1)+ (k12 + 2.0*(k22+k32) + k42)/6.0
end do 
end subroutine rk4_2d