program explicit
      parameter( nmax=1000 )
      parameter( mmax=100 )
      real u(nmax+1,mmax+1)
      real x(mmax+1), y(mmax+1)
      real alpha
*-- 出力ファイルをダンプ (保存) するための命令 ---
      CALL SWpSET( 'LDUMP', .TRUE. )
      CALL SWpSET( 'LWAIT', .FALSE. )
      CALL SWpSET( 'LWAIT1', .FALSE. )

*-- 定数の宣言 ---
      pi=3.1415926
*-- 空間領域の設定 ---
      xmin=0.0
      xmax=2.0*pi
*-- 空間きざみ幅の設定 ---
      dx=(xmax-xmin)/mmax
*-- 時間ステップの設定 ---
      dt=1.0/nmax
      alpha=(dt/dx)**2
*-- 位置座標の設定 ---
         do 230 j=1,mmax+1
            x(j)=xmin+dx*(j-1)
  230 continue

*-- 初期条件の設定 ---
      do 10 i=1,mmax+1
         u(1,i)=exp(-2.0*(x(i)-pi)**2)
  10  continue
      do 11 i=1,mmax+1
         u(2,i)=u(1,i)
  11  continue

*-- グラフの出力装置の指定 (IWS=4 は GTK 指定) ----
      WRITE(*,*) ' WORKSTATION ID (I)  ? ;'
      CALL SGPWSN
*      IWS=4
      READ (*,*) IWS

*-- 出力装置の呼び出し
      CALL GROPN( IWS )

*-- 時間ステップの計算 ---
      do 20 i=2,nmax
*-- 境界条件の設定 ---
         u(i+1,1)=0.0
         u(i+1,mmax+1)=0.0
         y(1)=u(i+1,1)
         y(mmax+1)=u(i+1,mmax+1)
         do 30 j=2,mmax
            u(i+1,j)=2.0*u(i,j)-u(i-1,j)
     &               +32.0*alpha*(u(i,j+1)-2.0*u(i,j)+u(i,j-1))
            y(j)=u(i+1,j)
  30  continue
*-- t=100 ごとで描画する設定 ---
      if(mod((i-1),50).eq.0)then
         CALL GRFRM
         CALL USSPNT( mmax+2, x, y )
         CALL GRSWND( xmin, xmax, -1.0, 1.0 )
         CALL GRSVPT(  0.2,  0.8,  0.2,  0.8 )
         CALL USPFIT
         CALL GRSTRF
         CALL USSTTL( 'X-AXIS', '', 'Y-AXIS', '' )
         CALL USDAXS
         CALL UULIN( mmax+1, x, y )
      end if
  20  continue

      stop
      end