      SUBROUTINE SIMP(FUN,A,B,C,N,RES)   

C-----THIS SUBROUTINE COMPUTES INTEGRALS OF USER SUPPLIED 
C-----FUNCTIONS BY SIMPSON'S RULE.  THE INPUTS ARE: 
C-----FUN = THE FUNCTION TO BE INTEGRATED 
C-----A = LOWER LIMIT 
C-----B = UPPER LIMIT 
C-----C = PARAMETER NEEDED TO DEFINE FUNCTION (COULD 
C-----    BE CHANGED TO AN ARRAY IF NECESSARY.) 
C-----N = NUMBER OF PANELS TO BE USED (MUST BE EVEN) 
C-----THE OUTPUT IS RES, THE RESULT. 

      EXTERNAL FUN 

C-----SET UP PARAMETERS 
      DEL = (B - A) / FLOAT(N) 
      M = N / 2	 
C-----FIRST COMPUTE THE ENDPOINT VALUES 
      RES = DEL * (FUN(A,C) + FUN(B,C))    
C-----NOW ADD THE ODD POINTS 
      DO 1 K=1,M 
      X = A + (2. * K - 1) * DEL 
      RES = RES + 4. * DEL * FUN(X,C) 
    1 CONTINUE  
C-----NOW ADD THE EVEN POINTS 
      MM1 = M - 1 
      DO 2 K=1,MM1 
      X = A + (2. * K) * DEL 
      RES = RES + 2. * DEL * FUN(X,C) 
    2 CONTINUE 
      RES = RES / 3. 

      END
