How to solve the following 3-dimensional system of non-linear inequalities in R?

I need to solve the following system of non-linear inequalities:

x*y >= 420
x*z >= 14
x*y*z < 5000

I have tried to find a similar question/solution or package that helps, but I struggle to apply it to that specific case.

The expected outcome should be a list of tuples (x,y,z). In the end, a 3-dimensional plot would be awesome but not really necessary (though should be easy to do it as soon as the solution list exists).

Edit: x, y and z are positive integers.

NB: It was recently clarified that the solutions should consist of integers. I'll give a solution for integers first, then the original answer without that condition.

The last condition x*y*z < 5000 indicates that all of the values must be less than 5000. That makes an exhaustive search feasible.

The way I'd do it is to set y to all possible values (i.e. 1:4999), then for each value find all possible x values. That will give a long list of (x,y) pairs. Then for each z value, select the pairs that meet all the conditions. For example:

x0 <- 1:4999
y0 <- 1:4999
z0 <- 1:4999

xy <- matrix(numeric(), ncol = 2, nrow = 0)
for (y in y0) {
  keep <- x0*y >= 420 & x0*y < 5000
  if (any(keep))
    xy <- rbind(xy, cbind(x0[keep], y))
}

xyz <- matrix(numeric(), ncol = 3, nrow = 0)
for (z in z0) {
  keep <- xy[, 1]*z >= 14 & xy[,1]*xy[,2]*z < 5000
  if (any(keep))
    xyz <- rbind(xyz, cbind(xy[keep, , drop=FALSE], z))
}

str(xyz)
#>  num [1:61833, 1:3] 420 421 422 423 424 425 426 427 428 429 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : NULL
#>   ..$ : chr [1:3] "" "y" "z"
head(xyz)
#>          y z
#> [1,] 420 1 1
#> [2,] 421 1 1
#> [3,] 422 1 1
#> [4,] 423 1 1
#> [5,] 424 1 1
#> [6,] 425 1 1
tail(xyz)
#>              y  z
#> [61828,] 2 222 11
#> [61829,] 2 223 11
#> [61830,] 2 224 11
#> [61831,] 2 225 11
#> [61832,] 2 226 11
#> [61833,] 2 227 11

^{Created on 2022-12-09 with reprex v2.0.2}

This is the original answer, without restricting to integer values:

The system of inequalities you give has an infinite number of solutions, so it's not possible to list them all in an R list. However, you could construct functions that describe the region holding the solutions.

One way to do that is to create two functions, with the first giving the range of y consistent with a particular x, and the second giving the range of z consistent with an (x, y) pair. For example,

min_y <- function(x) {
  420/x
}

range_z <- function(x, y) {
  lower <- 14/x
  upper <- 5000/x/y
  if (lower <= upper)
    c(lower, upper)
  else
    NULL
}

x <- 100
min_y(x)
#> [1] 4.2
range_z(x, min_y(x))
#> [1]  0.14000 11.90476
range_z(x, min_y(x) + 1)
#> [1] 0.140000 9.615385
range_z(x, min_y(x) + 1000)
#> NULL

^{Created on 2022-12-08 with reprex v2.0.2}

This shows that if x is 100, y must be bigger than 4.2. If you set y to 4.2 with x = 100, then z can be anything from 0.14 to 11.90476. If you set y to 5.2, z has to be between 0.14 and 9.615385. And if you set y to 1004.2, there are no solutions for z at all.