      PROGRAM BDYVAL 
C-----THIS PROGRAM SOLVES THE PROBLEM OF NEUTRON TRANSPORT IN A 
C-----FUEL ELEMENT.  SEE TEXT FOR DETAILS. 
      DIMENSION THETA(51),A(51),B(51),C(51),R(51),X(51),EXACT(51) 
C-----REQUEST INPUT 
      WRITE (6,100) 
  100 FORMAT ( ' TYPE ALPHA') 
      READ (5,110) ALPHA 
  110 FORMAT (F) 
      WRITE (6,120) 
  120 FORMAT ( ' TYPE BETA') 
      READ (5,110) BETA 
   10 WRITE (6,130) 
  130 FORMAT ( ' TYPE THE NUMBER OF INTERVALS') 
      READ (5,140) N 
  140 FORMAT (I) 
      NP1 = N + 1 
      DX = 1./FLOAT(N) 
C-----SET UP THE MATRIX FOR THE INTERIOR POINTS 
      DO 1 J=1,NP1 
      A(J) = 1. 
      B(J) = -2. - ALPHA*ALPHA*DX*DX 
      C(J) = 1. 
      R(J) = -DX*DX 
      X(J) = (J-1) * DX - .5 
      EXACT(J) = -((BETA * COSH(ALPHA*X(J))) / (ALPHA * SINH(.5*ALPHA) 
     1 + BETA * COSH(.5*ALPHA)) - 1.) / (ALPHA * ALPHA) 
    1 CONTINUE 
C-----NOW SET UP THE BOUNDARY CONDITIONS 
      C(1) = A(1) + C(1) 
      B(1) = B(1) - 2.* BETA * DX * A(1) 
      A(NP1) = A(NP1) + C(NP1) 
      B(NP1) = B(NP1) - 2.* BETA * DX * C(NP1) 
C-----SOLVE THE PROBLEM USING THE TRIDIAGONAL SYSTEM SOLVER; SEE 
C-----APPENDIX B FOR DETAILS. 
      CALL TRDIAG (NP1,A,B,C,THETA,R) 
C-----COMPUTE THE ERROR AND WRITE THE RESULTS 
      WRITE (6,150)  
  150 FORMAT ( 5H X(J),9X,5HEXACT,9X,5HTHETA,9X,5HERROR) 
      DO 2 J=1,NP1 
      ERROR = EXACT(J) - THETA(J) 
      WRITE (6,160) X(J),EXACT(J),THETA(J),ERROR 
  160 FORMAT ( 1F6.2,3X,3(E12.6,3X)) 
    2 CONTINUE 
      WRITE (6,170) 
  170 FORMAT ( ' ANOTHER CASE?, IF SO, TYPE 1')
      READ (5,140) IFF
      IF (IFF.EQ.1) GO TO 10
      STOP
      END
