C September 15, 1998
C ChE 515 illustration of solution of a linear BVP. 
C Driver program for LAPACK routine dgesv to solve A. x = b.
C
C Equation solved: d^2 (y)/dx^2 - (phi)^2 y = 0
C                  y(x=+-1) = 0
C Exact solution:  y = cosh(phi y)/cosh (phi).
C
C The problem represents the reaction diffusion problem 
C for a flat catalyst. 
C Discretization technique: Finite Differences.
C
      implicit double precision (a-h,o-z)
      parameter (neq=11, neqint=neq-2) 
C
C neq = total number of equations resulting from discretization.
C neqint = number of internal nodes (i.e., excluding the boundary nodes).
C
      dimension coef(neq,neq), rhs(neq) !Coefficient matrix and rhs vector.
      integer ipvt(neq) !Required for solver to store pivoting information.
C
C Specify discretization     
C
      delta = 2.d0/dfloat(neq-1)
C
C Thiele modulus
C 
      phi = 1.d0 
C
C Construct coefficient matrix and rhs vector. 
C
      do i=1,neq
         rhs(i) = 0.d0
        do j=1,neq
         coef(i,j) = 0.d0
        enddo
      enddo
c
      rhs(1) = 1.d0
      rhs(neq) = 1.d0
      coef(1,1) = 1.d0
      coef(neq,neq) = 1.d0
c
      do i=2,neq-1 
       coef(i,i-1) = 1.d0
       coef(i,i)   = -2.d0-phi*phi*delta*delta
       coef(i,i+1) = 1.d0  
      enddo
c
C Call dgesv to solve the system of equations
c
C ! !!!!! Use the tridiagonal system solver for problem 6 !!!!!!!
c 
      call dgesv(neq,1,coef,neq,ipvt,rhs,neq,info)
c
C Print results and compare with exact solution
c
      denom = (dexp(phi)+dexp(-phi))/2.d0
      do i=1,neq
       y = -1.d0+(i-1)*delta
       write(30,*) y, rhs(i)
       write(31,*) y, (dexp(phi*y)+dexp(-phi*y))/2.d0/denom
      enddo
C
      stop
      end
C
C Description of the routine dgesv. 
c
C     SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
*
*  -- LAPACK driver routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     March 31, 1993
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDB, N, NRHS
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
*     ..
*
*  Purpose
*  =======
*
*  DGESV computes the solution to a real system of linear equations
*     A * X = B,
*  where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
*
*  The LU decomposition with partial pivoting and row interchanges is
*  used to factor A as
*     A = P * L * U,
*  where P is a permutation matrix, L is unit lower triangular, and U is
*  upper triangular.  The factored form of A is then used to solve the
*  system of equations A * X = B.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The number of linear equations, i.e., the order of the
*          matrix A.  N >= 0.
*
*  NRHS    (input) INTEGER
*          The number of right hand sides, i.e., the number of columns
*          of the matrix B.  NRHS >= 0.
*
*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
*          On entry, the N-by-N coefficient matrix A.
*          On exit, the factors L and U from the factorization
*          A = P*L*U; the unit diagonal elements of L are not stored.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  IPIV    (output) INTEGER array, dimension (N)
*          The pivot indices that define the permutation matrix P;
*          row i of the matrix was interchanged with row IPIV(i).
*
*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
*          On entry, the N-by-NRHS matrix of right hand side matrix B.
*          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B.  LDB >= max(1,N).
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = i, U(i,i) is exactly zero.  The factorization
*                has been completed, but the factor U is exactly
*                singular, so the solution could not be computed.
*