!*************************************************************
!*   Differential equation y"=f(x,y,y') by Stormer's method  *
!* --------------------------------------------------------- *
!* SAMPLE RUN:                                               *
!* Integrate y" = 8yy / (1 + 2x) from x=0 to x=1,            *
!* with initial conditions: x(0)=0, y(0)=1 and y'(0)=-2      *
!* and compare with exact solution: y = 1 / (1 + 2x)         *
!*                                                           *
!* Output file (stormer.lst):                                *
!*                                                           *
!* --------------------------------------------------------- *
!*   Differential equation y"=f(x,y,y') by Stormer's method  *       
!* --------------------------------------------------------- *
!*        X           Y           Y exact         Error      *
!*     0.000       1.000000       1.000000     0.0000000000  *
!*     0.010       0.980392       0.980392     0.0000000001  *
!*     0.020       0.961538       0.961538     0.0000000295  *
!*     0.030       0.943396       0.943396     0.0000000457  *
!*     0.040       0.925926       0.925926     0.0000000974  *
!*     0.050       0.909091       0.909091     0.0000001285  *
!*     ...         ...            ...          ...           *
!*     0.950       0.344866       0.344828     0.0000381695  *
!*     0.960       0.342505       0.342466     0.0000388874  *
!*     0.970       0.340176       0.340136     0.0000396196  *
!*     0.980       0.337878       0.337838     0.0000403406  *
!*     0.990       0.335612       0.335570     0.0000410721  *
!*     1.000       0.333375       0.333333     0.0000418231  *
!*                                                           *
!*  End of file.                                             *
!*                        F90 Release By J-P Moreau, Paris.  *
!*                               (www.jpmoreau.fr)           *
!*************************************************************
  Program stormer
 
  real*8 c(4),x(4),y(4),z(4)
  real*8 h, a1, a2, a3
  real*8 F, Fx, G

  h=0.01d0 
  a1=1.08333333333333d0
  a2=-2.d0*(a1-1.d0)
  a3=a1-1.d0
  
  open(2,file='stormer.lst',status='unknown')
  Write(2,*) '-----------------------------------------------------------------'
  Write(2,*) '  Differential equation y"=f(x,y,y'') by Stormer''s method       '
  Write(2,*) '-----------------------------------------------------------------'
  !initial conditions
  x(1)=0.d0; y(1)=1.d0; z(1)=-2.d0
  yex=Fx(x(1))
  Write(2,*)'       X           Y           Y exact         Error'
  Write(2,20) x(1),y(1),yex,er
  do k=1,2
    ! call Runge-Kutta for first 2 steps
    call RK4D2(x(k),y(k),z(k),h,x(k+1),y(k+1),z(k+1))
    yex=Fx(x(k+1)); er=abs(yex-y(k+1));
    Write(2,20) x(k+1),y(k+1),yex,er
  end do
! main Stormer loop
10 continue
  do k=2,4 
    c(k)=G(x(5-k),y(5-k),z(5-k))
  end do
  y(4)=2.d0*y(3)-y(2)+h*h*(a1*c(2)+a2*c(3)+a3*c(4))
  x(4)=x(3)+h; yex=Fx(x(4)); er=abs(yex-y(4))
  write(2,20) x(4),y(4),yex,er
  do k=1,3
    x(k)=x(k+1); y(k)=y(k+1)
  end do
  if (x(3)<1.d0) goto 10  !end x value = 1.d0
  Write(2,*) ' '
  Write(2,*) ' End of file.'
  Close(2)
  print *,' '
  Write(*,*) ' Results in file stormer.lst'
  print *,' '
 
  stop

20 format('    ',f6.3,2f15.6,'  ',f15.10)

  END

     real*8 function F(x,y,z)
       real*8 x,y,z
       F=z
       return
     end

     real*8 function G(x,y,z)
       real*8 x,y,z
       G=8.d0*y*y/(1.d0+2.d0*x)
       return
     end

!    exact solution 
     real*8 function Fx(x)
       real*8 x
       Fx=1.d0/(1.d0+2.d0*x)
       return
     end

     SUBROUTINE RK4D2(x,y,z,h,x1,y1,z1)
       real*8 x,y,z,h,x1,y1,z1
       real*8 c1,d1,c2,d2,c3,d3,c4,d4
       real*8 F,G
       c1=F(x,y,z)
       d1=G(x,y,z)
       c2=F(x+h/2,y+h/2*c1,z+h/2*d1)
       d2=G(x+h/2,y+h/2*c1,z+h/2*d1)
       c3=F(x+h/2,y+h/2*c2,z+h/2*d2)
       d3=G(x+h/2,y+h/2*c2,z+h/2*d2)
       c4=F(x+h,y+h*c3,z+h*d3)
       d4=G(x+h,y+h*c3,z+h*d3)
       x1=x+h
       y1=y+h*(c1+2*c2+2*c3+c4)/6.
       z1=z+h*(d1+2*d2+2*d3+d4)/6.
       return
     end

!End of file stormer.f90