2階微分方程式の初等的な数値解法

!
!  y'' + p y' + qy = r
!    y2 = y
!    y1 = y2'
!  y1' + p y1 + qy2 = r
!  y1           = y2'
!
!	y1' = r - p y1 - q y2 = f1
!	y2' = y1            = f2
!
!
!
implicit none

 real :: dt
 integer :: n

 integer :: i
 real    :: dx, x, y1, y2, xmax, xmin
 real    :: f1, f2

 n = 1000		! number of dividing
 xmin = 0
 xmax = 10*3.1415926535
 dx = (xmax - xmin)/n

! initial value
 y1 = 0
 y2 = 0

 do i=1, n
  x = xmin + dx*i
  y1 = y1 + dx * f1(x,y1,y2)
  y2 = y2 + dx * f2(x,y1,y2)

  write(*,*) x, y2

 end do

end


real function f1(x, y1, y2)
implicit none
 real :: x, y1, y2

  f1 = sin(x) - 2*y1 - 2*y2
  return
end


real function f2(x, y1, y2)
implicit none
 real :: x, y1, y2

  f2 = y1
  return
end