Why is my approximation too large using Composite Simpson's rule in R (numerical integration)?

I am trying to approximate the following integral, using numerical integration in R:

,

where the function mu is defined by this formula:

To do this, I have implemented the Composite Simpson's rule as a function in R, which takes as parameters a function (integrand), the integration interval ([a,b]) and the number of subintervals desired (n).

I have tested my code on various different mathematical functions, and it seems to be working just fine. However, when I try to approximate the integral shown in the picture, the approximation becomes to large.

My method has been to first define the inner integral in terms of its Composite Simpson approximation as a function of t in R. Then, use the Composite Simpson's rule again, in order to calculate the outer integral by viewing the inner approximation as the function to be integrated.

When doing this, the inner approximation is correct when calculated by itself, as expected, but the approximation of the entire expression becomes too large, and I can't seem to figure out why.

I am comparing the approximations to those given by Maple; the inner expression calculated by itself, using t=20, should give 0.8157191, and the entire expression should be 12.837. R correctly calculates 0.8157191, but gives 32.9285 for the entire expression.

I have tried simplifying using numerous different mathematical functions, and making the functions independent of t in R, but all seems to result in the same error. So, to sum things up, my question is, why is only the outer integral being approximated wrongly?

I would be greatly appreciative of any hints or pointers - I have included my code illustrating the problem here:

compositesimpson <- function(integrand, a, b, n) {
  
  h<- (b-a)/n #THE DEFINITE INTERVAL IS SCALED BY
  #THE DESIRED NUMBER OF SUBINTERVALS
  
  xi<- seq.int(a, b, length.out = n+1) #DIVIDES THE DEFINITE INTERVAL INTO THE
  xi<- xi[-1]                          #DESIRED NUMBER OF SUBINTERVALS, 
  xi<- xi[-length(xi)]                 #EXCLUDING a AND b
  
  #THE APPROXIMATION ITSELF
  approks<- (h/3)*(integrand(a) + 2*sum(integrand(xi[seq.int(2, length(xi), 2)])) + 
                     4*sum(integrand(xi[seq.int(1, length(xi), 2)])) + integrand(b))
  return(approks)
  
}

# SHOULD YIELD -826.5755 BY Maple, SO THE FUNCTION IS WORKING HERE
ftest<- function(x) {
  return(exp(2*x)*sin(3*x))
}
compositesimpson(ftest, -4, 4, 100000)


# MU FUNCTION FOR TESTING
mu.01.kvinde<- function(x){ 0.000500 + 10^(5.728 + 0.038*(x+48) -10)}


#INNER INTEGRAL AS A FUNCTION OF ITS APPROXIMATION
indreintegrale.person1<- function(t){
  indre<- exp(-compositesimpson(mu.01.kvinde, 0, t, 100000))
  return(indre)
}

indreintegrale.person1(20) #YIELDS 0.8157191, WHICH IS CORRECT

compositesimpson(indreintegrale.person1, 20, 72, 100000) #YIELDS 32.9285,
#BUT SHOULD BE 12.837 ACCORDING TO MAPLE

This is something to do with trying to use vectorisation at two levels of recursion and it's not doing what you want it to. E.g. compare

indreintegrale.person1(20) 
#> [1] 0.8157191
indreintegrale.person1(c(20, 72))
#> [1] 0.8157191 0.4801160
indreintegrale.person1(72) 
#> [1] 2.336346e-10

I think the middle answer is wrong, but the other two are right.

Quickest fix, make this replacement:

indreintegrale.person1 <- function(t){
  sapply(t, function(t2) exp(-compositesimpson(mu.01.kvinde, 0, t2, 100000)))
}

and it now gives the answer you expect (but takes a bit longer to calculate!).