Issues with reproducing equation using the integral function in R

I am reading the following paper: https://www.uio.no/studier/emner/matnat/math/STK9190/v18/hjort-petrone-final.pdf and am trying to reproduce Figure 1 (currently just the left side). The left side of plot used equation 2 on page 3, i.e.,

which I have coded up in R as follows:

# Define the function to be integrated
integrand <- function(u, a, x) {
  denominator <- gamma(a - a * x) * gamma(a * x)
  numerator <- gamma(a)
  exponent <- a - a * x - 1
  result <- (numerator / denominator) * u^exponent * (1 - u)^(a * x - 1)
  return(result)
}

# Define the function representing the integral
integral <- function(x, y, a) {
  result <- integrate(integrand, lower = 0, upper = 1 - y, a = a, x = x)
  return(result$value)
}

# Define the partial derivative function
partial_derivative <- function(x, y, a, h = 1e-6) {
  (integral(x + h, y, a) - integral(x, y, a)) / h
}

# Test the partial derivative function
x <- 0
y <- 0.75
a <- 0.1
result <- partial_derivative(x, y, a)
print(result)

However, this seems to always result in an error when varying the values of a and x. The errors tend to vary which is why I am not posting it all, but here is an example:

Error in integrate(integrand, lower = 0, upper = y - 1, a = a, x = x) : 
  non-finite function value

I am hoping that it is something simple that somewhat can spot right away as I can't seem to find it. Likewise, I am open to other solutions for calculating equation (2) if there is a better solution.

This is a beta distribution. Compare the following with your function:

integral_x <- function(x, y, a) pbeta(1 - y, a * (1 - x), a*x)

integral_x(x=0.1, y=0.2, a=0.3)
[1] 0.1288844
integral(x=0.1, y=0.2, a=0.3) # from your function above
[1] 0.1288844

The function gives exact results as your function. You can then carry out numerical differentiation on the integral_x function to obtain the needed derivative:

partial_dev <- function(x, y, a, h = 1e-4){
  (integral_x(x+h, y, a) - integral_x(x, y, a))/h
}

partial_dev(x = 0, y = 0.75, a = 0.3)
[1] 0.7043913