c From csnet!mit-multics.arpa!UBC.mailnet!USER=NBAF Thu, 5 Feb 87 21:54:24 PST
C  PROBLEM 1 - SEE COMPANION PAPER [1] 
C 
                IMPLICIT REAL*8 (A-H,O-Z) 
      REAL*8 FSPACE(2000), ZETA(4), TOL(2), Z(4), U(4), ERR(4) 
      INTEGER ISPACE(200), M(1), IPAR(11), LTOL(2) 
      EXTERNAL FSUB, DFSUB, GSUB, DGSUB, DUMMY 
      REAL*8 FIXPNT(1)
C 
      WRITE (6,99) 
   99 FORMAT(1H1, 35H EXAMPLE OF A SIMPLE PROBLEM SETUP. 
     .         /  46H  UNIFORMLY LOADED BEAM OF VARIABLE STIFFNESS, 
     .         /  32H  SIMPLY SUPPORTED AT BOTH ENDS. /) 
C 
C     ONE DIFFERENTIAL EQUATION OF ORDER 4. 
      M(1) = 4 
C     GIVE LOCATION OF BOUNDARY CONDITIONS 
      ZETA(1) = 1.D0 
      ZETA(2) = 1.D0 
      ZETA(3) = 2.D0 
      ZETA(4) = 2.D0 
C     SET UP PARAMETER ARRAY. 
C     USE DEFAULT VALUES FOR ALL PARAMETERS EXCEPT FOR INITIAL 
C     MESH SIZE, NO. OF TOLERANCES AND SIZES OF WORK ARRAYS 
      DO 10 I=1,11 
   10   IPAR(I) = 0 
      IPAR(3) = 1 
      IPAR(4) = 2 
      IPAR(5) = 2000 
      IPAR(6) = 200 
      IPAR(7) = 1
C     TWO ERROR TOLERANCES (ON U AND ITS SECOND DERIVATIVE) 
      LTOL(1) = 1 
      LTOL(2) = 3 
      TOL(1) = 1.D-11 
      TOL(2) = 1.D-11 
C 
      CALL TIME(0)
      CALL COLSYS (1, M, 1.D0, 2.D0, ZETA, IPAR, LTOL, TOL, 
     .             FIXPNT, ISPACE, FSPACE, IFLAG, FSUB, 
     .             DFSUB, GSUB, DGSUB, DUMMY) 
      CALL TIME(3,1)
C 
      IF (IFLAG .NE. 1)  STOP 
      CALL TIME(0)
C     CALCULATE THE ERROR AT 101 POINTS USING THE KNOWN 
C     EXACT SOLUTION 
      X = 1.D0 
      DO 20 I=1,4 
   20    ERR(I) = 0.D0 
      DO 40 J=1,101 
         CALL APPSLN (X, Z, FSPACE, ISPACE) 
         CALL EXACT (X, U) 
         DO 30 I=1,4 
   30       ERR(I) = DMAX1(ERR(I), DABS(U(I)-Z(I))) 
   40    X = X + .01D0 
      WRITE(6,100) (ERR(I),I=1,4) 
  100 FORMAT(/27H ERROR TOLERANCES SATISFIED//22H THE EXACT ERRORS ARE, 
     .       / 7X,4D12.4) 
      CALL TIME(3,1)
      STOP 
      END 
                SUBROUTINE FSUB (X, Z, F) 
      REAL*8 Z(4), F(1), X 
      F(1) = (1.D0 - 6.D0*X**2*Z(4) - 6.D0*X*Z(3)) / X**3 
      RETURN 
      END 
                SUBROUTINE DFSUB (X, Z, DF) 
      REAL*8 Z(4), DF(1,4), X 
      DF(1,1) = 0.D0 
      DF(1,2) = 0.D0 
      DF(1,3) = -6.D0/X**2 
      DF(1,4) = -6.D0/X 
      RETURN 
      END 
                SUBROUTINE GSUB (I, Z, G) 
      REAL*8 Z(4), G 
      GO TO (1, 2, 1, 2), I 
    1 G = Z(1) - 0.D0 
      RETURN 
    2 G = Z(3) - 0.D0 
      RETURN 
      END 
                SUBROUTINE DGSUB (I, Z, DG) 
      REAL*8 Z(4), DG(4) 
      DO 10 J=1,4 
   10    DG(J) = 0.D0 
      GO TO (1, 2, 1, 2), I 
    1 DG(1) = 1.D0 
      RETURN 
    2 DG(3) = 1.D0 
      RETURN 
      END