Actually to understand the possible answers for the given question we can try to solve the system of inequalities graphically.

There was a nice answer concerning plotting of inequations in R at stackowerflow. Using the given aproach we can plot the following graph:

library(ggplot2)

fun1 <- function(x) 2*x - 3        # this is the same as -2x + y <= -3
fun2 <- function(x) -1.25*x + 2.5  # 1.25x + y <= 2.5
fun3 <- function(x) -3             # y >= -3
x1 = seq(-1,5, by = 1/16)
mydf = data.frame(x1, y1=fun1(x1), y2=fun2(x1),y3= fun3(x1))
mydf <-  transform(mydf, z = pmax(y3,pmin(y1,y2)))
ggplot(mydf, aes(x = x1)) + 
  geom_line(aes(y = y1), colour = 'blue') +
  geom_line(aes(y = y2), colour = 'green') +
  geom_line(aes(y = y3), colour = 'red') +
  geom_ribbon(aes(ymin=y3,ymax = z), fill = 'gray60')

All the possible (infinite by number) tuples lie inside the gray triangle. The vertexes can be found using the following code.

obj <- numeric(2)
mat <- matrix(c(-2, 1.25, 1, 1), nrow = 2)
rhs <- matrix(c(-3, 2.5), nrow = 2)

aPoint <- solve(mat, rhs)

mat <- matrix(c(-2, 0, 1, 1), nrow = 2)
rhs <- matrix(c(-3, -3), nrow = 2)

bPoint <- solve(mat, rhs)

mat <- matrix(c(1.25, 0, 1, 1), nrow = 2)
rhs <- matrix(c(2.5, -3), nrow = 2)

cPoint <- solve(mat, rhs) 

Note the order of arguments of matrices.

And you get the coordinates:

> aPoint
          [,1]
[1,] 1.6923077
[2,] 0.3846154
> bPoint
     [,1]
[1,]    0
[2,]   -3
> cPoint
     [,1]
[1,]  4.4
[2,] -3.0

All the codes below are with base R only (no need library(Rglpk))

1. Corner Points

If you want to get all the corner points, here is one option

A <- matrix(c(-2, 1.25, 0, 1, 1, -1), nrow = 3)

b <- c(-3, 2.5, 3)

# we use `det` to check if the coefficient matrix is singular. If so, we return `Inf`.
xh <-
  combn(nrow(A), 2, function(k) {
    if (det(A[k, ]) == 0) {
      rep(NA, length(k))
    } else {
      solve(A[k, ], b[k])
    }
  })

# We filter out the points that satisfy the constraint
corner_points <- t(xh[, colSums(A %*% xh <= b, na.rm = TRUE) == length(b)])

such that

> corner_points
         [,1]       [,2]
[1,] 1.692308  0.3846154
[2,] 0.000000 -3.0000000
[3,] 4.400000 -3.0000000

2. Possible Tuples

If you want to have multiple tuples, e.g., n=10, we can use Monte Carlo simulation (based on the obtained corner_points in the previous step) to select the tuples under the constraints:

xrange <- range(corner_points[, 1])
yrange <- range(corner_points[, 2])
n <- 10
res <- list()
while (length(res) < n) {
  px <- runif(1, xrange[1], xrange[2])
  py <- runif(1, yrange[1], yrange[2])
  if (all(A %*% c(px, py) <= b)) {
    res[length(res) + 1] <- list(c(px, py))
  }
}

and you will see n possible tuples in a list like below

> res
[[1]]
[1]  3.643167 -2.425809

[[2]]
[1]  2.039007 -2.174171

[[3]]
[1]  0.4990635 -2.3363637

[[4]]
[1]  0.6168402 -2.6736421

[[5]]
[1]  3.687389 -2.661733

[[6]]
[1]  3.852258 -2.704395

[[7]]
[1] 1.7571062 0.1067597

[[8]]
[1]  3.668024 -2.771307

[[9]]
[1]  2.108187 -1.365349

[[10]]
[1]  2.106528 -2.134310

First of all, the matrix representing the three equations needs a small correction, because R fills matrices column by column :

  -2x + y <= -3
1.25x + y <= 2.5
        y >= -3

mat <- matrix(c(-2, 1.25, 0, 1, 1, 1), nrow = 3
# and not : mat <- matrix(c(-2, 1, 1.25, 1, 0, 1), nrow = 3)

To get different tuples, you could modify the objective function :

obj <- numeric(2) results in an objective function 0 * x + 0 * y which is always equal to 0 and can't be maximized : the first valid x,y will be selected.

Optimization on x is achieved by using obj <- c(1,0), resulting in maximization / minimization of 1 * x + 0 * y.
Optimization on y is achieved by using obj <- c(0,1).

#setting the bounds is necessary, otherwise optimization occurs only for x>=0 and y>=0
bounds <- list(lower = list(ind = c(1L, 2L), val = c(-Inf, -Inf)),
               upper = list(ind = c(1L, 2L), val = c(Inf, Inf)))

# finding maximum x: obj = c(1,0), max = T
Rglpk_solve_LP(obj = c(10,0), mat = mat, dir = dir, rhs = rhs,bound=bounds, max = T)$solution
# [1] 4.4 -3.0

# finding minimum x: obj = c(1,0), max = F
Rglpk_solve_LP(obj = c(10,0), mat = mat, dir = dir, rhs = rhs,bound=bounds, max = F)$solution
#[1]  0 -3

# finding maximum y: obj = c(0,1), max = T
Rglpk_solve_LP(obj = c(0,1), mat = mat, dir = dir, rhs = rhs,bound=bounds, max = T)$solution
#[1] 1.6923077 0.3846154