C 下の例では q = 0, g = 1 となっている．
C 区間 [0,1] 上で -u''+qu = g を斉次ディリクレ条件の下で解く
      PROGRAM SOLVEQ
      PARAMETER (NMAX=9999)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(NMAX,NMAX),F(NMAX),X(NMAX)
      WRITE(*,'(1H ,A5,A5,A23,A18)')
     +        'N','H','Numerical solution','Exact solution'
      N=10
      DO 100 J=1,4
        H=1.0D0/N
        CALL SETMAT(A,H,N-1)
        CALL SETVEC(F,H,N-1)
        CALL SOLVE(A,F,X,N-1)
        WRITE(*,'(I6,F8.5,2F20.15)')N,H,X(N/2),SOLTN(DBLE(1)/2)
        N=N*10
  100 CONTINUE
      END
C
C ガウスの消去法
C
      SUBROUTINE SOLVE(A,F,X,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(N,N),F(N),X(N)
C 前進消去
      DO 300 I=1,N
         IMAX=I
         DO 50 K=I+1,N
            IF (DABS(A(IMAX,I)).LT.DABS(A(K,I)))IMAX=K !ピボットを探す
   50    CONTINUE
         IF (IMAX.GT.I) CALL SWAP(A,F,I,IMAX,N)  ! ピボットを現在行に移す
         DO 200 J=I+1,N
            FACTOR=A(J,I)/A(I,I)
            DO 100 K=I+1,N
                A(J,K)=A(J,K)-A(I,K)*FACTOR
  100       CONTINUE
            F(J)=F(J)-F(I)*FACTOR
  200    CONTINUE
C         DO 250 K=I+1,N           ! 実際には使わないので不要
C            A(K,I)=0.0D0
C  250    CONTINUE
  300 CONTINUE
C 後退代入
      DO 600 I=N,1,-1
          X(I)=F(I)
          DO 500 J=I+1,N
              X(I)=X(I)-X(J)*A(I,J)
  500     CONTINUE
          X(I)=X(I)/A(I,I)
  600 CONTINUE
      RETURN
      END
C
      SUBROUTINE SWAP(A,F,I,J,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(N,N),F(N)
C      WRITE(*,*)'SWAP',I,' AND ',J
      DO 100 K=I,N
          DUMMY=A(I,K)
          A(I,K)=A(J,K)
          A(J,K)=DUMMY
  100 CONTINUE
      DUMMY=F(I)
      F(I)=F(J)
      F(J)=DUMMY
      RETURN
      END
C
C 微分方程式の離散化から行列を作る
C
      SUBROUTINE SETMAT(A,H,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(N,N)
      DO 100 I=1,N
         DO 100 J=1,N
            A(I,J)=0.0D0
  100 CONTINUE
      X=H
      A(1,1)=2/H/H+Q(X)
      X=X+H
      DO 200 I=2,N
         A(I,I)=2/H/H+Q(X)
         A(I-1,I)=-1/H/H
         A(I,I-1)=-1/H/H
         X=X+H
  200 CONTINUE
      RETURN
      END
C
C 右辺の既知函数から離散化ベクトルを作る
C
      SUBROUTINE SETVEC(F,H,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION F(N)
      X=H
      DO 100 I=1,N
          F(I)=G(X)
          X=X+H
  100 CONTINUE
      RETURN
      END
C
C 右辺の函数の定義
C
      DOUBLE PRECISION FUNCTION G(X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      G=1.0D0
      G=-(2-12*X+12*X**2)
      RETURN
      END
C
C 係数の函数の定義
C
      DOUBLE PRECISION FUNCTION Q(X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      Q=0.0D0
      RETURN
      END

C 真の解の定義
C
      DOUBLE PRECISION FUNCTION SOLTN(X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SOLTN=X*(1-X)/2
      SOLTN=X*X*(1-X)*(1-X)
      RETURN
      END