Q:

The code runs successfully but no graphs as an output

I am trying to replicate the work of Nonlinear-Infection by Fabian Dablander in the case of Pakistan. But when I run the codes below then it doesn't make any graph as an output"
plot_SIR <- function(res, main = '') {
cols <- brewer.pal(3, 'Set1')
matplot(
res, type = 'l', col = cols, axes = FALSE, lty = 1, lwd = 2,
ylab = 'Subpopulations(t)', xlab = 'Time t', xlim = c(0, 4000),
ylim = c(0, 1), main = main, cex.main = 1.75, cex.lab = 1.5,
font.main = 1, xaxs = 'i', yaxs = 'i'
)

axis(1, cex.axis = 1.25)
axis(2, las = 2, cex.axis = 1.25)
legend(
3000, 0.65, col = cols, legend = c('S', 'I', 'R'),
lty = 1, lwd = 2, bty = 'n', cex = 1.5
)
}

A:

Taking a quick look, you need to run the other function to get the param: res, then stick it inside the original function just as you described in the post:
solve_SIR <- function(S0, I0, beta = 1, gamma = 1, delta_t = 0.01, times = 8000) {
  res <- matrix(
    NA, nrow = times, ncol = 4, dimnames = list(NULL, c('S', 'I', 'R', 'Time'))
  )
  res[1, ] <- c(S0, I0, 1 - S0 - I0, delta_t)

  dS <- function(S, I) -beta * I * S
  dI <- function(S, I)  beta * I * S - gamma * I

  for (i in seq(2, times)) {
    S <- res[i-1, 1]
    I <- res[i-1, 2]

    res[i, 1] <- res[i-1, 1] + delta_t * dS(S, I)
    res[i, 2] <- res[i-1, 2] + delta_t * dI(S, I)
    res[i, 4] <- delta_t * i
  }

  res[, 3] <- 1 - res[, 1] - res[, 2]
  res
}

use function to gen param input and choosing random values for S0 and I0 this will output the graph you desired.
plot_SIR(solve_SIR(0.95, 0.05))

