C  PCPRB7.FOR  27 February 1991
C
      PROGRAM PCPRB7
C
C  Materially nonlinear problem
C
C  Reference:
C
C  Ivo Babuska, Werner C Rheinboldt,
C  Reliable Error Estimations and Mesh adaptation for the Finite Element Method,
C  In: Computation Methods in Nonlinear Mechanics,
C  Edited by J T Oden,
C  North Holland Publishing Company,
C  Amsterdam 1980, Pages 67-108.
C
C  The function:
C
C  The general problem is the two point boundary value problem:
C
C  -D/DX PHI(X,U,DU/DX,T)  +  PSI(X,U,DU/DX,T) = 0,
C
C  U(X=0)=U0, U(X=1)=U1.
C
C  The particular problem solved here has the form:
C
C  U'' + T*SIN(U+U**2+U**3)=0
C
C  U(0)=U(1)=0.0
C
C  U is approximated by dividing the interval into NINT subintervals,
C  and using the NPOLYS Legendre functions of orders 0 through NPOLYS-1
C  on each subinterval.  The NINT*NPOLYS coefficients of these basis
C  functions form the representation of U.  
C
C  Moreover, there are NINT*NPOLYS finite element equations set up,
C  of the form 
C
C  INTEGRAL(Subinterval I) -PHI(X,U,U',T)*PL'(J) + PSI(X,U,U',T)*PL(J) = 0
C
C  There is a complication however, in that Lagrange multipliers are
C  used to enforce continuity of the solution U at the interface
C  nodes between neighboring subintervals.  These conditions create
C  new variables, and new equations, and even modify the above standard
C  finite element equations.
C
C  Options:
C
C  ICHOOZ=1  Piecewise linear, 1 continuity condition
C  ICHOOZ=2  Piecewise cubic, 1 continuity condition
C  ICHOOZ=3  Piecewise cubic, 2 continuity conditions
C  ICHOOZ=4  Piecewise quintic, 1 continuity condition
C  ICHOOZ=5  Piecewise quintic, 2 continuity conditions
C  ICHOOZ=6  Piecewise quintic, 3 continuity conditions
C
C  All options use 8 intervals.
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   MINT
      INTEGER   MVAR
C
      PARAMETER (MINT=8)
      PARAMETER (MVAR=MINT*6+2+(MINT-1)*3)
      PARAMETER (LIW=MVAR+29)
      PARAMETER (LRW=29+(6+MVAR)*MVAR)
C
      EXTERNAL  BVAL
      EXTERNAL  DENSLV
      EXTERNAL  FP0011
      EXTERNAL  FX0011
      EXTERNAL  PITCON
      EXTERNAL  UVAL
C
      INTRINSIC MOD
C
      REAL      BCONE
      REAL      BCZERO
      REAL      DBCODT
      REAL      DBCZDT
      REAL      FPAR(1)
      REAL      GCOEF
      REAL      GPOINT
      INTEGER   I
      INTEGER   ICHOOZ
      INTEGER   IERROR
      INTEGER   IHOLD
      INTEGER   II
      INTEGER   IPAR(1)
      INTEGER   ISKIP
      INTEGER   IWORK(LIW)
      INTEGER   LOUNIT
      CHARACTER NAME*12
      INTEGER   NBCZ
      INTEGER   NBCO
      INTEGER   NDRV
      INTEGER   NINT
      INTEGER   NPOLYS
      INTEGER   NVAR
      INTEGER   NVARY
      INTEGER   NVARZ
      REAL      PL
      REAL      PLD
      REAL      RWORK(LRW)
      REAL      TEMP
      REAL      THETA
      REAL      U
      REAL      UPRYM
      REAL      XG
      REAL      XL
      REAL      XR(MVAR)
      REAL      XRR
C
      COMMON /AUXMEM/ BCZERO(8),DBCZDT(8),BCONE(8),DBCODT(8),
     1 NPOLYS,NDRV,NVARY,NVARZ,NBCZ,NBCO,
     2 PL(8),PLD(8),GCOEF(8),GPOINT(8),NINT,THETA(10,10)
C
      DATA GCOEF /
     1  .1012285363,  .2223810345,  .3137066459,  .3626837834,
     2  .3626837834,  .3137066459,  .2223810345,  .1012285363/
      DATA GPOINT /
     1  .9602898565,  .7966664774,  .5255324099,  .1834346425,
     2 -.1834346425, -.5255324099, -.7966664774, -.9602898565/
      DATA THETA/1.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0, 0.0,
     3 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     5 1.0, 3.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     7 1.0, 6.0, 15.0, 15.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     9 1.0, 10.0, 45.0, 105.0, 105.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     2 1.0, 15.0, 105.0, 420.0, 945.0, 945.0, 0.0, 0.0, 0.0, 0.0,
     4 1., 21., 210., 1260., 4725., 10395., 10395., 0., 0., 0.,
     6 1., 28., 378., 3150., 17325., 62370., 135135., 135135., 0., 0.,
     8 1., 36., 630., 6930., 51975., 
     * 270270., 945945., 2027025., 2027025., 0.,
     1 1., 45., 990., 13860., 135135., 
     * 945945., 4729725., 16216200.,34459425.,34459425./
C
      LOUNIT=6
C
C  Zero out work arrays
C
      DO 10 I=1,LIW
        IWORK(I)=0
10      CONTINUE
      DO 20 I=1,LRW
        RWORK(I)=0.0
20      CONTINUE
C
C  Set number of variables
C
      ICHOOZ=1
      NINT=8
      IF(ICHOOZ.EQ.1)NPOLYS=2
      IF(ICHOOZ.EQ.2.OR.ICHOOZ.EQ.3)NPOLYS=4
      IF(ICHOOZ.EQ.4.OR.ICHOOZ.EQ.5.OR.ICHOOZ.EQ.6)NPOLYS=6
      NVARY=NINT*NPOLYS
      NBCZ=1
      IF(ICHOOZ.EQ.1.OR.ICHOOZ.EQ.2.OR.ICHOOZ.EQ.4)NDRV=1
      IF(ICHOOZ.EQ.3.OR.ICHOOZ.EQ.5)NDRV=2
      IF(ICHOOZ.EQ.6)NDRV=3
      NBCO=1
      NVARZ=NBCZ+(NINT-1)*NDRV+NBCO
      NVAR=NVARY+NVARZ+1
C
C  Starting point
C
      DO 110 I=1,NVAR
        XR(I)=0.0   
110     CONTINUE
C
C  We perturb the zero solution, by setting the 5-th coefficient
C  to some nonzero value.  We will then also require that the program
C  find an exact solution by correcting this starting point, but 
C  without altering the 5-th coefficient, thus guaranteeing that the
C  solution is not all zero.
C
      IHOLD=5
      XR(IHOLD)=0.001
C
C  IWORK(1)=0     ; This is a startup
C  IWORK(2)=IHOLD ; Use index IHOLD for first parameter
C  IWORK(3)=0     ; Program may choose index
C  IWORK(4)=0     ; Update jacobian every newton step
C  IWORK(5)=0     ; Seek no target values.
C  IWORK(6)=NVAR  ; Seek limit points in index NVAR
C  IWORK(7)=1     ; Moderate amount of output
C  IWORK(8)=6     ; Output unit to use
C  IWORK(9)=0     ; Use user's jacobian routine
C
C  IWORK(17)=20   ; (Unusual setting)  Increase the allowed number of Newton
C                   iterations to 20.
C
      IWORK(1)=0
      IWORK(2)=IHOLD
      IWORK(3)=0
      IWORK(4)=0
      IWORK(5)=0
      IWORK(6)=NVAR
      IWORK(7)=1
      IWORK(8)=LOUNIT
      IWORK(9)=0
      IWORK(17)=20
C
C  RWORK(1)=0.0001 ; Absolute error tolerance
C  RWORK(2)=0.0001 ; Relative error tolerance
C  RWORK(3)=0.5    ; Minimum stepsize
C  RWORK(4)=0.5    ; Maximum stepsize
C  RWORK(5)=0.25   ; Starting stepsize
C  RWORK(6)=1.0    ; Starting direction
C  RWORK(7)=0.0    ; Target value 
C
      RWORK(1)=0.0001
      RWORK(2)=0.0001
      RWORK(3)=0.001
      RWORK(4)=0.5
      RWORK(5)=0.25
      RWORK(6)=1.0
      RWORK(7)=0.0
C
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'PCPRB7'
      WRITE(LOUNIT,*)'PITCON sample program'
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'The materially nonlinear problem.'
      WRITE(LOUNIT,*)' '
      IF(ICHOOZ.EQ.1)THEN
        WRITE(LOUNIT,*)'Piecewise linears'
      ELSEIF(ICHOOZ.EQ.2.OR.ICHOOZ.EQ.3)THEN
        WRITE(LOUNIT,*)'Piecewise cubics'
      ELSE
        WRITE(LOUNIT,*)'Piecewise quintics.'
        ENDIF
      IF(ICHOOZ.EQ.6)THEN
        WRITE(LOUNIT,*)'3 continuity conditions at breakpoints'
      ELSEIF(ICHOOZ.EQ.5.OR.ICHOOZ.EQ.3)THEN
        WRITE(LOUNIT,*)'2 continuity conditions at breakpoints'
      ELSE
        WRITE(LOUNIT,*)'1 continuity condition at breakpoints.'
        ENDIF
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Step  Type of Point     Lambda'
      WRITE(LOUNIT,*)' '
      I=0
      NAME='Start point  '
      WRITE(LOUNIT,'(1X,I3,2X,A12,2X,3G14.6)')I,NAME,XR(NVAR)
C
C  Take another step.
C
      DO 30 I=1,30
C
        CALL PITCON(FP0011,FPAR,FX0011,IERROR,IPAR,IWORK,LIW,
     *  NVAR,RWORK,LRW,XR,DENSLV)
C
        IF(IWORK(1).EQ.1)THEN
          NAME='Corrected    '
        ELSEIF(IWORK(1).EQ.2)THEN
          NAME='Continuation '
        ELSEIF(IWORK(1).EQ.3)THEN
          NAME='Target point '
        ELSEIF(IWORK(1).EQ.4)THEN
          NAME='Limit point  '
          ENDIF
C
        WRITE(LOUNIT,'(1X,I3,2X,A12,2X,3G14.6)')I,NAME,XR(NVAR)
C
C  Print out shape of rod every 5 steps
C
        IF(MOD(I,5).EQ.0)THEN
          WRITE(LOUNIT,*)' '
          WRITE(LOUNIT,*)'Current rod coordinates:'
          WRITE(LOUNIT,*)' '
          WRITE(LOUNIT,*)'X,  U(X),  U''(X)'
          WRITE(LOUNIT,*)' '
          DO 25 II=1,NINT
            TEMP=NINT
            XL=(II-1)/TEMP
            XRR=(II)/TEMP
            ISKIP=(II-1)*NPOLYS
            XG=0.5*(XL+XRR)
            CALL BVAL(NPOLYS,PL,PLD,XG,XL,XRR)
            CALL UVAL(NPOLYS,PL,PLD,XG,U,UPRYM,XR,ISKIP)
            WRITE(LOUNIT,'(1X,3G14.6)')XG,U,UPRYM
25          CONTINUE
          WRITE(LOUNIT,*)' '
          ENDIF
C
        IF(IWORK(1).EQ.3)GO TO 40
C
        IF(IERROR.NE.0)THEN
          WRITE(LOUNIT,*)'Iteration halted, IERROR=',IERROR
          GO TO 40
          ENDIF
C
30      CONTINUE
40    CONTINUE
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Jacobians      ',IWORK(19)
      WRITE(LOUNIT,*)'Factorizations ',IWORK(20)
      WRITE(LOUNIT,*)'Solves         ',IWORK(21)
      WRITE(LOUNIT,*)'Functions      ',IWORK(22)
      STOP
      END
      SUBROUTINE FX0011(NVAR,FPAR,IPAR,X,FX,IERROR)
C
C***********************************************************************
C
      EXTERNAL  BD0011
      EXTERNAL  BVAL
      EXTERNAL  FVAL
      EXTERNAL  UVAL
C
      INTEGER   NVAR
C
      REAL      BCONE
      REAL      BCZERO
      REAL      DBCODT
      REAL      DBCZDT
      REAL      DTDX
      REAL      DTDXL
      REAL      DTDXR
      REAL      FPAR(*)
      REAL      FX(NVAR)
      REAL      GCOEF
      REAL      GPOINT
      REAL      H2I
      REAL      H2IL
      REAL      H2IR
      INTEGER   I
      INTEGER   IEQN
      INTEGER   IERROR
      INTEGER   INDX
      INTEGER   IPAR(*)
      INTEGER   ISKIP
      INTEGER   IVAR
      INTEGER   J
      INTEGER   K
      INTEGER   KHI
      INTEGER   L
      INTEGER   LHIL
      INTEGER   LHIR
      INTEGER   NBCZ
      INTEGER   NBCO
      INTEGER   NCL
      INTEGER   NCR
      INTEGER   NDRV
      INTEGER   NDSUM
      INTEGER   NEQN
      INTEGER   NINT
      INTEGER   NODES
      INTEGER   NPOLYS
      INTEGER   NPSUM
      INTEGER   NVARY
      INTEGER   NVARZ
      REAL      PHI
      REAL      PHIPT
      REAL      PHIPU
      REAL      PHIPUP
      REAL      PL
      REAL      PLD
      REAL      PSI
      REAL      PSIPT
      REAL      PSIPU
      REAL      PSIPUP
      REAL      S
      REAL      TCON
      REAL      TEMP
      REAL      TERM
      REAL      TERML
      REAL      TERMR
      REAL      THETA
      REAL      U
      REAL      UPRYM
      REAL      X(NVAR)
      REAL      XC
      REAL      XG
      REAL      XL
      REAL      XR
C
      COMMON /AUXMEM/ BCZERO(8),DBCZDT(8),BCONE(8),DBCODT(8),
     1 NPOLYS,NDRV,NVARY,NVARZ,NBCZ,NBCO,
     2 PL(8),PLD(8),GCOEF(8),GPOINT(8),NINT,THETA(10,10)
C
C  Zero out FX
C
      TCON=X(NVAR)
      CALL BD0011
      NEQN=NVAR-1
      DO 10 J=1,NEQN
   10   FX(J)=0.0
C
C  1. SET UP THE TERMS (A*Y) INVOLVING THE BIVARIATE FORM
C
C  I  FOR EACH INTERVAL
C
      DO 40 I=1,NINT
        TEMP=NINT
        XL=(I-1)/TEMP
        XR=(I)/TEMP
        ISKIP=(I-1)*NPOLYS
C
C  J  FOR EACH GAUSS POINT, EVALUATE THE INTEGRAND
C
        DO 30 J=1,8
          XG=XL+(GPOINT(J)+1.0)*(XR-XL)/2.0
          CALL BVAL(NPOLYS,PL,PLD,XG,XL,XR)
          CALL UVAL(NPOLYS,PL,PLD,XG,U,UPRYM,X,ISKIP)
          CALL FVAL(XG,U,UPRYM,TCON,PHI,PHIPU,PHIPUP,PHIPT,
     *              PSI,PSIPU,PSIPUP,PSIPT)
C
C  L   PROJECT ONTO EACH TEST FUNCTION PL(L) AND PLD(L)
C
          IEQN=ISKIP
          DO 20 L=1,NPOLYS
            IEQN=IEQN+1
            TERM=GCOEF(J)*(XR-XL)*(PSI*PL(L)+PHI*PLD(L))/2.0
            FX(IEQN)=FX(IEQN)+TERM
   20       CONTINUE
   30     CONTINUE
   40   CONTINUE
C
C  2. ADD THE TERMS INVOLVING THE CONTINUITY OF THE TEST FUNCTIONS
C  WHICH ARE THE TERMS B*Z IN  F=A*Y + B*Z
C
      NCL=NVARY
      NCR=NVARY+NBCZ
C
C  I  FOR EACH INTERVAL
C
      DO 110 I=1,NINT
        TEMP=NINT
        XL=(I-1)/TEMP
        XR=(I)/TEMP
        DTDX=2.0/(XR-XL)
C
C  J   FOR THE POLYNOMIALS USED IN APPROXIMATING EACH U,
C  COUNT CONDITIONS AT LEFT ENDPOINT, LHIL, AND AT RIGHT, LHIR
C  IF WE ARE IN THE FIRST OR LAST INTERVAL, ONE OF
C  THESE WILL BE BOUNDARY CONDITIONS
C
        LHIL=NDRV
        LHIR=NDRV
        IF (I.EQ.1) LHIL=NBCZ
        IF (I.EQ.NINT) LHIR=NBCO
        IF (LHIL.LE.0.AND.LHIR.LE.0) GO TO 100
C
C  K  FOR EACH TEST FUNCTION PL(K)
C
        DO 90 K=1,NPOLYS
          TERML=0.0
          TERMR=0.0
          S=(-1.0)**(K+1)
          IEQN=(I-1)*NPOLYS+K
C
C  L  FOR DERIVATIVES 0 THRU NBCZ SPECIFIED AT 0.0 TO BE ZERO
C  OR DERIVATIVES 0 THRU NDRV SPECIFIED TO BE CONTINUOUS
C  AT INTERIOR NODES
C  OR DERIVATIVES 0 THRU NBCO SPECIFIED AT 1.0 TO BE ZERO
C
          IF (LHIL.LE.0) GO TO 60
          H2I=1.0
          DO 50 L=1,LHIL
            S=-S
            INDX=NCL+L
            TERML=TERML+S*X(INDX)*H2I*THETA(L,K)
            H2I=H2I*DTDX
   50       CONTINUE
   60     IF (LHIR.LE.0) GO TO 80
          H2I=1.0
          DO 70 L=1,LHIR
            INDX=NCR+L
            TERMR=TERMR+H2I*THETA(L,K)*X(INDX)
            H2I=H2I*DTDX
   70       CONTINUE
   80     FX(IEQN)=FX(IEQN)+TERML+TERMR
   90     CONTINUE
        IF (LHIL.GT.0) NCL=NCL+LHIL
        IF (LHIR.GT.0) NCR=NCR+LHIR
  100   CONTINUE
  110   CONTINUE
C
C  3. CREATE THE TERMS FOR THE U FUNCTIONS AND THEIR DERIVATIVES
C  THE MATRIX TERMS  ( C*Y )
C  ONE EQUATION IS GENERATED FOR COMPONENT AND CONDITION
C
      NPSUM=0
      DTDXR=0.0
      DTDXL=0.0
C
C  I  FOR EACH NODE
C
      NDSUM=NVARY
      NODES=NINT+1
      DO 210 I=1,NODES
        TEMP=NINT
        IF (I.GT.1) XL=(I-2)/TEMP
        XC=(I-1)/TEMP
        IF (I.LT.NODES) XR=(I)/TEMP
        IF (XC.NE.XL) DTDXL=2.0/(XC-XL)
        IF (XR.NE.XC) DTDXR=2.0/(XR-XC)
        H2IL=1.0
        H2IR=1.0
C
C  K   FOR DERIVATIVES 0 THRU NBCZ SPECIFIED AT 0.0
C  OR DERIVATIVES 0 THRU NDRV SPECIFIED TO BE CONTINUOUS
C  AT INTERIOR NODES
C  OR DERIVATIVES 0 THRU NBCO SPECIFIED AT 1.0
C
        KHI=NDRV
        IF (I.EQ.1) KHI=NBCZ
        IF (I.EQ.NODES) KHI=NBCO
        IF (KHI.LE.0) GO TO 190
        DO 180 K=1,KHI
          S=(-1.0)**(K+1)
C
C  L   SET UP THE TERM FROM THE LEFT HAND INTERVAL
C
          TERML=0.0
          IF (I.GT.1) GO TO 120
          TERML=BCZERO(K)
          GO TO 140
  120     DO 130 L=1,NPOLYS
            IVAR=NPSUM+L-NPOLYS
            TERML=TERML+X(IVAR)*H2IL*THETA(K,L)
  130       CONTINUE
C
C  L   SET UP THE TERM FROM THE RIGHT HAND INTERVAL
C
  140     IF (I.LT.NODES) GO TO 150
          TERMR=-BCONE(K)
          GO TO 170
  150     TERMR=0.0
          DO 160 L=1,NPOLYS
            IVAR=NPSUM+L
            S=-S
            TERMR=TERMR+S*X(IVAR)*H2IR*THETA(K,L)
  160       CONTINUE
  170     IEQN=NDSUM+K
          FX(IEQN)=TERML+TERMR
          H2IL=H2IL*DTDXL
          H2IR=H2IR*DTDXR
  180     CONTINUE
        NDSUM=NDSUM+KHI
  190   NPSUM=NPSUM+NPOLYS
  200   CONTINUE
  210   CONTINUE
      RETURN
      END
      SUBROUTINE FP0011(NVAR,FPAR,IPAR,X,FPRIME,IERROR)
C
C***********************************************************************
C
C  JACOBIAN OF MATERIALLY NONLINEAR PROBLEM.
C
      EXTERNAL  BD0011
      EXTERNAL  BVAL
      EXTERNAL  FVAL
      EXTERNAL  UVAL
C
      INTEGER   NVAR
C
      REAL      BCONE
      REAL      BCZERO
      REAL      DBCODT
      REAL      DBCZDT
      REAL      DTDX
      REAL      FPAR(*)
      REAL      FPRIME(NVAR,NVAR)
      REAL      GCOEF
      REAL      GPOINT
      REAL      H2I
      INTEGER   I
      INTEGER   IEQN
      INTEGER   IERROR
      INTEGER   IPAR(*)
      INTEGER   ISKIP
      INTEGER   IVAR
      INTEGER   J
      INTEGER   K
      INTEGER   KHIL
      INTEGER   KHIR
      INTEGER   L
      INTEGER   LHIL
      INTEGER   LHIR
      INTEGER   N
      INTEGER   NBCZ
      INTEGER   NBCO
      INTEGER   NCL
      INTEGER   NCR
      INTEGER   NDRV
      INTEGER   NINT
      INTEGER   NPOLYS
      INTEGER   NPSUM
      INTEGER   NVARY
      INTEGER   NVARZ
      REAL      PHI
      REAL      PHIPT
      REAL      PHIPU
      REAL      PHIPUP
      REAL      PL
      REAL      PLD
      REAL      PSI
      REAL      PSIPT
      REAL      PSIPU
      REAL      PSIPUP
      REAL      S
      REAL      SUM
      REAL      TCON
      REAL      TEMP
      REAL      TERM
      REAL      TERM1
      REAL      TERM2
      REAL      TERM3
      REAL      TERM4
      REAL      THETA
      REAL      U
      REAL      UPRYM
      REAL      X(NVAR)
      REAL      XG
      REAL      XL
      REAL      XR
C
      COMMON /AUXMEM/ BCZERO(8),DBCZDT(8),BCONE(8),DBCODT(8),
     1 NPOLYS,NDRV,NVARY,NVARZ,NBCZ,NBCO,
     2 PL(8),PLD(8),GCOEF(8),GPOINT(8),NINT,THETA(10,10)
C
C  ZERO OUT THE MATRIX
C
      TCON=X(NVAR)
      CALL BD0011
      DO 20 I=1,NVAR
        DO 10 J=1,NVAR
          FPRIME(I,J)=0.0
   10     CONTINUE
   20   CONTINUE
C
C  1.  SET UP THE TERMS FROM THE BIVARIATE FORM (A*Y )
C
C      I  FOR EACH INTERVAL (XL,XR)
C
      DO 80 I=1,NINT
        TEMP=NINT
        XL=(I-1)/TEMP
        XR=(I)/TEMP
        ISKIP=(I-1)*NPOLYS
C
C  J  FOR EACH GAUSS POINT XG IN THE INTERVAL
C  NOTE THAT THE LEGENDRE POLYNOMIALS AND DERIVATIVES ARE SET UP
C  BY THE CALL TO VALUE
C
        DO 70 J=1,8
          XG=XL+0.5*(XR-XL)*(GPOINT(J)+1.0)
          CALL BVAL(NPOLYS,PL,PLD,XG,XL,XR)
          CALL UVAL(NPOLYS,PL,PLD,XG,U,UPRYM,X,ISKIP)
          CALL FVAL(XG,U,UPRYM,TCON,PHI,PHIPU,PHIPUP,PHIPT,
     *              PSI,PSIPU,PSIPUP,PSIPT)
C
C  L  FOR EACH LEGENDRE POLYNOMIAL COEFFICIENT
C
          IEQN=ISKIP
          DO 50 L=1,NPOLYS
            IEQN=IEQN+1
            TERM=GCOEF(J)*0.5*(XR-XL)*(PSIPT*PL(L)+PHIPT*PLD(L))
            FPRIME(IEQN,NVAR)=FPRIME(IEQN,NVAR)+TERM
C
C  N   FOR EACH Y-COEFFICIENT OF A U
C
            DO 30 N=1,NPOLYS
              IVAR=NPOLYS*(I-1)+N
              TERM1=PSIPU*PL(N)*PL(L)
              TERM2=PSIPUP*PLD(N)*PL(L)
              TERM3=PHIPU*PL(N)*PLD(L)
              TERM4=PHIPUP*PLD(N)*PLD(L)
              SUM=TERM1+TERM2+TERM3+TERM4
              FPRIME(IEQN,IVAR)=FPRIME(IEQN,IVAR)+
     *        GCOEF(J)*0.5*(XR-XL)*SUM
   30         CONTINUE
   50       CONTINUE
   70     CONTINUE
   80   CONTINUE
C
C  2. ADD THE TERMS INVOLVING THE CONTINUITY OF THE TEST FUNCTIONS
C  WHICH ARE THE TERMS B*Z IN  F=A*Y+B*Z
C
C
C     I  FOR EACH INTERVAL,
C
      DO 150 I=1,NINT
        NCL=NVARY
        IF (I.GT.1) NCL=NVARY+NBCZ+(I-2)*NDRV
        NCR=NVARY+NBCZ+(I-1)*NDRV
        TEMP=NINT
        XL=(I-1)/TEMP
        XR=(I)/TEMP
        DTDX=2.0/(XR-XL)
C
C  J   FOR THE POLYNOMIALS USED IN APPROXIMATING EACH U,
C      COUNT CONDITIONS AT LEFT ENDPOINT, LHIL, AND AT RIGHT, LHIR
C      IF WE ARE IN THE FIRST OR LAST INTERVAL, ONE OF
C      THESE WILL BE BOUNDARY CONDITIONS
C
          LHIL=NDRV
          LHIR=NDRV
          IF (I.EQ.1) LHIL=NBCZ
          IF (I.EQ.NINT) LHIR=NBCO
          IF (LHIL.LE.0.AND.LHIR.LE.0) GO TO 130
C
C  K  FOR EACH TEST FUNCTION PL(K)
C
          DO 120 K=1,NPOLYS
            S=(-1.0)**(K+1)
            IEQN=(I-1)*NPOLYS+K
C
C  L  FOR DERIVATIVES 0 THRU NBCZ SPECIFIED AT 0.0 TO BE ZERO
C     OR DERIVATIVES 0 THRU NDRV SPECIFIED TO BE CONTINUOUS
C     AT INTERIOR NODES
C     OR DERIVATIVES 0 THRU NBCO SPECIFIED AT 1.0 TO BE ZERO
C
C     EVALUATE CONTRIBUTION FROM LEFT ENDPOINT
C
            IF (LHIL.LE.0) GO TO 100
            H2I=1.0
            DO 90 L=1,LHIL
              S=-S
              IVAR=NCL+L
              FPRIME(IEQN,IVAR)=S*H2I*THETA(L,K)
              H2I=H2I*DTDX
   90         CONTINUE
C
C  EVALUATE CONTRIBUTION FROM RIGHT ENDPOINT
C
  100       IF (LHIR.LE.0) GO TO 120
            H2I=1.0
            DO 110 L=1,LHIR
              IVAR=NCR+L
              FPRIME(IEQN,IVAR)=H2I*THETA(L,K)
              H2I=H2I*DTDX
  110         CONTINUE
  120       CONTINUE
          IF (LHIL.GT.0) NCL=NCL+LHIL
          IF (LHIR.GT.0) NCR=NCR+LHIR
  130     CONTINUE
  140     CONTINUE
  150   CONTINUE
C
C  3. CREATE THE TERMS FOR THE U FUNCTIONS AND THEIR DERIVATIVES
C  THE MATRIX TERMS  ( C*Y )
C  ONE EQUATION IS GENERATED FOR COMPONENT AND CONDITION
C
C
C  I  FOR EACH INTERVAL
C
      DO 230 I=1,NINT
        NCL=NVARY
        IF (I.GT.1) NCL=NVARY+NBCZ+(I-2)*NDRV
        NCR=NVARY+NBCZ+(I-1)*NDRV
        NPSUM=(I-1)*NPOLYS
        TEMP=NINT
        XL=(I-1)/TEMP
        XR=(I)/TEMP
        DTDX=2.0/(XR-XL)
        H2I=1.0
C
C  K   FOR DERIVATIVES 0 THRU NBCZ SPECIFIED AT 0.0
C      OR DERIVATIVES 0 THRU NDRV SPECIFIED TO BE CONTINUOUS
C      AT INTERIOR NODES
C      OR DERIVATIVES 0 THRU NBCO SPECIFIED AT 1.0
C
        KHIL=NDRV
        IF (I.EQ.1) KHIL=NBCZ
        IF (KHIL.LE.0) GO TO 180
        DO 170 K=1,KHIL
          IEQN=NCL+K
C
C  L   SET UP THE TERM FROM THE LEFT HAND ENDPOINT
C
          IF (I.EQ.1) FPRIME(IEQN,NVAR)=DBCZDT(K)
          S=(-1.0)**(K+1)
          DO 160 L=1,NPOLYS
            IVAR=NPSUM+L
            S=-S
            FPRIME(IEQN,IVAR)=S*H2I*THETA(K,L)
  160       CONTINUE
          H2I=H2I*DTDX
  170     CONTINUE
        NCL=NCL+KHIL
  180   CONTINUE
        H2I=1.0
        KHIR=NDRV
        IF (I.EQ.NINT) KHIR=NBCO
        IF (KHIR.LE.0) GO TO 210
        DO 200 K=1,KHIR
          IEQN=NCR+K
          FPRIME(IEQN,NVAR)=0.0
          IF (I.EQ.NINT) FPRIME(IEQN,NVAR)=-DBCODT(K)
          DO 190 L=1,NPOLYS
            IVAR=NPSUM+L
            FPRIME(IEQN,IVAR)=H2I*THETA(K,L)
  190       CONTINUE
          H2I=H2I*DTDX
  200     CONTINUE
        NCR=NCR+KHIR
  210   NPSUM=NPSUM+NPOLYS
  230   CONTINUE
      RETURN
      END
      SUBROUTINE BVAL(NPOLYS,PL,PLD,XG,XL,XR)
C
C***********************************************************************
C
C  AUXILLIARY ROUTINE FOR MATERIALLY NONLINEAR PROBLEM.
C
C  EVALUATE LEGENDRE POLYNOMIALS AND DERIVATIVES
C
      REAL      A
      INTEGER   I
      INTEGER   NPOLYS
      REAL      PL(8)
      REAL      PLD(8)
      REAL      T
      REAL      XG
      REAL      XL
      REAL      XR
C
      T=-1.0+2.0*(XG-XL)/(XR-XL)
      PL(1)=1.0
      PLD(1)=0.0
      PL(2)=T
      PLD(2)=1.0
      A=0.0
      DO 10 I=3,8
        A=A+1.0
        PL(I)=((A+A+1.0)*T*PL(I-1)-A*PL(I-2))/(A+1.0)
        PLD(I)=((A+A+1.0)*(T*PLD(I-1)+PL(I-1))-A*PLD(I-2))/(A+1.0)
10      CONTINUE
C
      DO 30 I=1,8
        PLD(I)=2.0*PLD(I)/(XR-XL)
30      CONTINUE
      RETURN
      END
      SUBROUTINE UVAL(NPOLYS,PL,PLD,XG,U,UPRYM,X,ISKIP)
C
C  SUM ON COEFFICENTS TO EVALUATE U AND UPRYM
C
      INTEGER   NPOLYS
C
      INTEGER   INDX
      INTEGER   ISKIP
      INTEGER   J
      REAL      PL(NPOLYS)
      REAL      PLD(NPOLYS)
      REAL      U
      REAL      UPRYM
      REAL      X(*)
      REAL      XG
C
      U=0.0
      UPRYM=0.0
      DO 40 J=1,NPOLYS
        INDX=ISKIP+J
        U=U+X(INDX)*PL(J)
        UPRYM=UPRYM+X(INDX)*PLD(J)
40      CONTINUE
      RETURN
      END
      SUBROUTINE FVAL(XG,U,UPRYM,TCON,PHI,PHIPU,PHIPUP,PHIPT,
     * PSI,PSIPU,PSIPUP,PSIPT)
C
C***********************************************************************
C
C  AUXILLIARY ROUTINE FOR MATERIALLY NONLINEAR PROBLEM.
C
      INTRINSIC COS
      INTRINSIC SIN
C
      REAL      COS
      REAL      PHI
      REAL      PHIPU
      REAL      PHIPUP
      REAL      PHIPT
      REAL      PSI
      REAL      PSIPU
      REAL      PSIPUP
      REAL      PSIPT
      REAL      SIN
      REAL      TCON
      REAL      U
      REAL      UPRYM
      REAL      XG
C
      PHI=-UPRYM
      PHIPU=0.0
      PHIPUP=-1.0
      PHIPT=0.0
C
      PSI=TCON*SIN(U*(1.0+U*(1.0+U)))
      PSIPU=TCON*(1.0+U*(2.0+3.0*U))*COS(U*(1.0+U*(1.0+U)))
      PSIPUP=0.0
      PSIPT=SIN(U*(1.0+U*(1.0+U)))
C
      RETURN
      END
      SUBROUTINE BD0011
C
C***********************************************************************
C
C  AUXILLIARY ROUTINE SETS BOUNDARY VALUES
C  FOR MATERIALLY NONLINEAR PROBLEM.
C
      REAL      BCONE
      REAL      BCZERO
      REAL      DBCODT
      REAL      DBCZDT
      REAL      GCOEF
      REAL      GPOINT
      INTEGER   NBCZ
      INTEGER   NBCO
      INTEGER   NDRV
      INTEGER   NINT
      INTEGER   NPOLYS
      INTEGER   NVARY
      INTEGER   NVARZ
      REAL      PL
      REAL      PLD
      REAL      THETA
C
      COMMON /AUXMEM/ BCZERO(8),DBCZDT(8),BCONE(8),DBCODT(8),
     1 NPOLYS,NDRV,NVARY,NVARZ,NBCZ,NBCO,
     2 PL(8),PLD(8),GCOEF(8),GPOINT(8),NINT,THETA(10,10)
C
      BCONE(1)=0.0
      BCZERO(1)=0.0
      DBCODT(1)=0.0
      DBCZDT(1)=0.0
      RETURN
      END