      PROGRAM example                       ! 主プログラムの開始
*     初期値とパラメータの設定
      IMPLICIT NONE
      REAL*8 delta_x, x, xmax, y
      CALL initial(y,x,delta_x,xmax)
      DO WHILE(x .LE. xmax)
         CALL Euler(y,x,delta_x)            ! 簡単なオイラーのアルゴリズムを使用
         WRITE(*,'(2F13.6)') x, y
      END DO
      END
C
      SUBROUTINE initial(y,x,delta_x,xmax)
      IMPLICIT NONE
      REAL*8 y, x, delta_x, xmax
      x = 1                                 ! x の初期値
      y = 1                                 ! y の初期値
      xmax = 2                              ! x の最大値
      delta_x = 0.1D0                       ! 分割幅
      WRITE(*,'(2A13)') ' x ', ' y '        ! 見出し部分の表示
      WRITE(*,*)                            ! 空白行
      WRITE(*,'(2F13.6)') x, y              ! 初期値の表示
      END
C
      SUBROUTINE Euler(y,x,delta_x)
*     わかりやすさのためだけに変数 slope と change を導入する
      IMPLICIT NONE
      REAL*8 change, slope, y, x, delta_x
      slope = 2*x                           ! 区間の初めでの傾き
      change = slope*delta_x                ! 区間における変化の計算
      y = y + change                        ! y の新しい値
      x = x + delta_x                       ! 区間の終りでの x の値
      END
