program grank_nicolson
      parameter( nmax=10 )
      parameter( mmax=100 )
      real a, b, c
      real x(mmax+1), y(mmax+1), ca(mmax+1), da(mmax+1)
      real d(mmax+1)
      real u(nmax,mmax+1)
      real dt, dx, r
      real pi, xmin, xmax, kappa
      real tmin, tmax
      pi=3.14159            ! 円周率
      kappa=0.01             ! 拡散係数
*-- 出力ファイルをダンプ (保存) するための命令 ---
      CALL SWpSET( 'LDUMP', .TRUE. )
      CALL SWpSET( 'LWAIT', .FALSE. )
      CALL SWpSET( 'LWAIT1', .FALSE. )

*-- x 方向の描画領域 ---
      xmin=0.0
      xmax=1.0
*-- きざみ幅の設定 ---
      dx=(xmax-xmin)/mmax
      dt=0.1
      alpha=(kappa**2)*(dt/dx**2)
*-- 陰解法に使用する3重対角系行列の成分の設定 ---
      a=-alpha
      b=1.0+2.0*alpha
      c=a
*-- 連立方程式を計算するための係数の計算 ---
          ca(2)=c/b
        do 2 i=3,mmax
           ca(i)=c/(b-a*ca(i-1))
  2   continue
*-- 初期条件の設定 ---
      do 10 i=1,mmax+1
         x(i)=xmin+dx*(i-1)
         u(1,i)=exp(-20.0*(x(i)-0.5)**2)
  10  continue

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

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

*-- 時間発展の計算 ---
      do 30 i=1,nmax-1
*-- 境界条件の設定 ---
         u(i+1,1)=0.0
         u(i+1,mmax+1)=0.0
*-- 連立方程式を計算するための係数の計算その2 ---
         d(2)=u(i,2)+alpha*u(i+1,1)
         d(mmax)=u(i,mmax)+alpha*u(i+1,mmax+1)
         do 20 j=3,mmax-1
            d(j)=u(i,j)
  20  continue
         d(2)=d(2)/b
         do 21 j=3,mmax
            d(j)=d(j)/(b-a*ca(j-1))
  21  continue
*-- 後退代入の計算 ---
         u(i+1,mmax)=d(mmax)
         do 22 j=mmax-1,2,-1
            u(i+1,j)=d(j)-ca(j)*u(i+1,j+1)
  22  continue
         do 23 j=1,mmax+1
            y(j)=u(i+1,j)
  23  continue

         CALL GRFRM
         CALL USSPNT( mmax+1, x, y )
         CALL GRSWND( xmin, xmax, 0.0, 1.0 )
         CALL GRSVPT(  0.2,  0.8,  0.2,  0.8 )
         CALL USPFIT
         CALL GRSTRF
         CALL USSTTL( 'X-AXIS', '', 'TEMPERATURE', '' )
         CALL USDAXS
         CALL UULIN( mmax+1, x, y )

  30  continue
      stop
      end