付録2‐1ホインスキームのプログラム例
              program Heun
        !振動方程式のホインスキームの式
        !U* = U^(n) + iω Δ tU^(n)
        !U^(n + 1) = U^(n) + 1/2iω Δ t(U^(n) + U*)
        !初期値は1.0d0
        !振動数ω はπ 

!-------宣言文ここから--------
       implicit none

       real(8) omega                              !倍精度実数で振動数を宣言

       real(8) pi, h, time, amf, amfg             !倍精度実数でπ とΔ tと時間幅と増幅値を宣言

       real,allocatable,dimension(:) :: a, b      !1列配列を単精度で配列数は後で決めるように宣言 

       integer i, n, IWS                          !ループ整数と切断数とDCLの宣言

       complex(kind (0d0)) U_a, U_b, ei, kai !U_aはU^(n + 1),U_bはU^(n)でeiは虚数,kaiはU^*の宣言
!------宣言文ここまで---------

       pi = (4.0d0)*(atan(1.0d0)) !π の定義
       ei = (0.0d0, 1.0d0)        !虚数iの定義
       omega = pi                 !振動数の定義
       time = 2.0d0               !全体の時間

       write(*,*) "Enter in noumber"

       read(*,*) n                !切断数の入力

       allocate(a(n + 1))         !Re[U]の配列数をn + 1に定義
       allocate(b(n + 1))         !tの配列数をn + 1に定義

       h = time/dble(n)           !ステップ幅の定義

       open(10, file='Heunjissuu.dat', status='replace') !データの書き込みのための呼び出し

       U_b = (1.0d0, 0.0d0) !U^(n)の初期値
       kai = (1.0d0, 0.0d0) !U*の初期値
       U_a = (0.0d0, 0.0d0) !U^(n + 1)の初期値


       a(1) = real(U_b)  !Re[U]の配列の初期値

       b(1) = 0.0e0      !tの初期値
!------ホインスキーム初め----------------------------
       do i = 1, n
         kai = U_b + (h*ei*omega*U_b)                       !U*の計算
         U_a = U_b + 0.5d0*h*(ei*omega*kai + ei*omega*U_b)  !Uの計算

         amf = abs(1.0d0 + ei*omega*h - 0.5d0*(omega*h)**2.0d0) !増幅値の計算
         amfg = abs(U_a)/abs(U_b)                           !解析解の増幅値

         a(i + 1) = real(U_a)                               !Re[U]を配列に代入
         b(i + 1) = real(i)*real(h)                         !tを配列に代入

         write(10,*) real(U_a)                              !datへの書き込み
         write(*,*) amf, amfg                               !増幅値の書き出し
         write(*,*) U_a,  a(i + 1), b(i + 1), omega*h       !数値解と配列とpの書き出し

         U_b = U_a                                          !U^(n + 1)をU^(n)に代入

       end do
!-----ホインスキーム終わり------------------------
 
       close(10)
!------DCL初め--------------------------------------       
       CALL SGPWSN
       read(*,*) IWS

       call GROPN( IWS )
       call GRFRM
       
       call USSTTL( 't', 'x-unit', 'Re[U]', 'y-unit' )
       call USGRPH( n + 1, b, a )

       call GRCLS
!------DCL終わり------------------------------------
       stop

       end program Heun