Using antiD function for variance of gamma distribution

Question

This is my first post here and I hope I'll follow all the rules of the community.

I'm trying to calculate variance of gamma distribution with shape parameter 2 and scale parameter 3 in R using function antiD from mosaic package. The R code I use is the following

stopifnot(require(mosaic))

f <- function(y) {
     dgamma(y, shape = 2, scale = 3)
}

mean_integral <- antiD( z*f(z) ~ z )
mn <- mean_integral(10^4)

g <- function(y) {
     (y - mn)^2
}

variance <- antiD(f(x)*g(x) ~ x)
variance(10^5)
## [1] 7.115334e-09

The problem is that the number I get doesn't make sense as the variance for Gamma distribution with those parameters should be equal to 2*3^2 = 18 (Wiki page on Gamma distribution). Moreover if I put 10^4 as an upper bound (the default lower bound is 0) for variance() it will return the following:

variance(10^4)
## [1] 18

And the integral from 10^4 to 10^5 will be:

variance(10^5) - variance(10^4)
## [1] -18

Does anyone know why variance(10^5) produce nonsensical results in this case? I also will be grateful for any additional comments on the style of the post.