x'' = p(t)x'(t)+q(t)x(t)+r(t)
  x(b) = beta in [a,b]

program test_linsht
!
! Linear Shooting method
! To approximate the solution of boundary value problem
!           x'' = p(t)x'(t)+q(t)x(t)+r(t)
!'          x(b) = beta in [a,b]
! and using Runge-Kutta method of order N = 4
!
!
!
implicit none
integer,parameter :: n = 2
integer,parameter :: nstep = 20
real(8) :: U(nstep)             ! U contains the solutions at each step
real(8) :: V(nstep),out(nstep)  ! V just like U for another equation
real(8) :: a,b,h,alph,beta,w
real(8) :: x(0:n)
integer :: i
a = 0.d0
b = 4.0
alph = 125.d-2
beta  = -95.d-2
h = (b-a)/nstep
! x = t0 x0 x'0
x = (/0.d0, alph, 0.d0/)

call rks4U(n,h,x,nstep,U)
a = 0.d0
b = 4.0
x = (/0.d0, 0.d0, 1.d0 /)
h = (b-a)/nstep

call rks4V(n,h,x,nstep,V)

! calculate the solution to the boundary value broblem
w = (beta - U(nstep))/V(nstep)
write (*,*) w,U(nstep),V(nstep)
!w = 0.485884
out = U + w * V
do i = 1,nstep
    write(*,100) i,i*h,U(i), w*V(i), out(i)
enddo
100 format (1x,I5,f10.2,f15.6,f15.6,f15.6)
end program

subroutine xprsysU(n,x,f)
! x prime of system of equations

    real(8), dimension(0:n) :: x,f
    integer :: n
    ! time is introduced as a new differential equation
    f(0) = 1.d0
    f(1) = x(2)
    f(2) = 2.*x(0)/(1+x(0)*x(0))*x(2)- 2./(1.+x(0)*x(0))*x(1) + 1.0
end subroutine

subroutine xprsysV(n,x,f)
! x prime of system of equations

    real(8), dimension(0:n) :: x,f
    integer :: n
    ! time is introduced as a new differential equation
    f(0) = 1.d0
    f(1) = x(2)
    f(2) = 2.*x(0)/(1+x(0)*x(0))*x(2)- 2./(1.+x(0)*x(0))*x(1)
end subroutine



subroutine rks4U(n,h,x,nstep,xout)

implicit none
      real(8) :: x(0:n)
      real(8) :: y(0:n), f(0:n,4)
      real(8) :: xout(nstep)
      integer :: i,k,n,nstep
      real(8) :: h
      !write(*,200) "k","t","x","y"
      !write(*,100) 0,x
      f = 0.
out:  do k = 1,nstep
          call xprsysU(n,x,f(0,1))
in1:      do i = 0,n
              y(i) = x(i) + 0.5*h*f(i,1)
              !print *, f(i,1)
          end do in1
          call xprsysU(n,y,f(0:n,2))
in2:      do i = 0,n
              y(i) = x(i) + 0.5*h*f(i,2)
          end do in2
          call xprsysU(n,y,f(0:n,3))
in3:      do i = 0,n
              y(i) = x(i) + h*f(i,3)
          end do in3
          call xprsysU(n,y,f(0:n,4))
in4:      do i = 0,n
              x(i) = x(i) + (h/6.0)* (f(i,1) + 2.0*(f(i,2) + f(i,3)) + f(i,4))
          end do in4
          xout(k) = x(1)
      end do out
      !100 format (1x,I5,f10.2,f15.8,f15.8)
      !200 format (1x,A4,2A10,A16)
end subroutine rks4U


subroutine rks4V(n,h,x,nstep,xout)

implicit none
      real(8) :: x(0:n)
      real(8) :: y(0:n), f(0:n,4)
      real(8) :: xout(nstep)
      integer :: i,k,n,nstep
      real(8) :: h
      !write(*,200) "k","t","x","y"
      !write(*,100) 0,x
      f = 0.
out:  do k = 1,nstep
          call xprsysV(n,x,f(0,1))
in1:      do i = 0,n
              y(i) = x(i) + 0.5*h*f(i,1)
              !print *, f(i,1)
          end do in1
          call xprsysV(n,y,f(0:n,2))
in2:      do i = 0,n
              y(i) = x(i) + 0.5*h*f(i,2)
          end do in2
          call xprsysV(n,y,f(0:n,3))
in3:      do i = 0,n
              y(i) = x(i) + h*f(i,3)
          end do in3
          call xprsysV(n,y,f(0:n,4))
in4:      do i = 0,n
              x(i) = x(i) + (h/6.0)* (f(i,1) + 2.0*(f(i,2) + f(i,3)) + f(i,4))
          end do in4
          xout(k) = x(1)
      end do out
      !100 format (1x,I5,f10.2,f15.8,f15.8)
      !200 format (1x,A4,2A10,A16)
end subroutine rks4V

program test_linsht
   implicit none
   ...
   call rk4(n, h, x, nstep, xprsysU, U)
   call rk4(n, h, x, nstep, xprsysV, V)
   ...
contains

subroutine xprsysU(n,x,f)
! x prime of system of equations
   real(8), dimension(0:n) :: x,f
   integer :: n
   ! time is introduced as a new differential equation
   f(0) = 1.d0
   f(1) = x(2)
   f(2) = 2.*x(0)/(1+x(0)*x(0))*x(2)- 2./(1.+x(0)*x(0))*x(1) + 1.0
end subroutine

subroutine rk4(n, h, x, nstep, sub, xout)
   real(8) :: x(0:n)
   real(8) :: y(0:n), f(0:n,4)
   interface
      subroutine sub(n, x, f)
         real(8), dimension(0:n) :: x,f
         integer :: n
      end subroutine sub
   end interface
   ...
end subroutine
end program test_linsht