get x-value given y-value: general root finding for linear / non-linear interpolation function

First of all, let me copy in the stable solution for linear interpolation proposed in my previous answer.

## given (x, y) data, find x where the linear interpolation crosses y = y0
## the default value y0 = 0 implies root finding
## since linear interpolation is just a linear spline interpolation
## the function is named RootSpline1
RootSpline1 <- function (x, y, y0 = 0, verbose = TRUE) {
  if (is.unsorted(x)) {
     ind <- order(x)
     x <- x[ind]; y <- y[ind]
     }
  z <- y - y0
  ## which piecewise linear segment crosses zero?
  k <- which(z[-1] * z[-length(z)] <= 0)
  ## analytical root finding
  xr <- x[k] - z[k] * (x[k + 1] - x[k]) / (z[k + 1] - z[k])
  ## make a plot?
  if (verbose) {
    plot(x, y, "l"); abline(h = y0, lty = 2)
    points(xr, rep.int(y0, length(xr)))
    }
  ## return roots
  xr
  }

For cubic interpolation splines returned by stats::splinefun with methods "fmm", "natrual", "periodic" and "hyman", the following function provides a stable numerical solution.

RootSpline3 <- function (f, y0 = 0, verbose = TRUE) {
  ## extract piecewise construction info
  info <- environment(f)$z
  n_pieces <- info$n - 1L
  x <- info$x; y <- info$y
  b <- info$b; c <- info$c; d <- info$d
  ## list of roots on each piece
  xr <- vector("list", n_pieces)
  ## loop through pieces
  i <- 1L
  while (i <= n_pieces) {
    ## complex roots
    croots <- polyroot(c(y[i] - y0, b[i], c[i], d[i]))
    ## real roots (be careful when testing 0 for floating point numbers)
    rroots <- Re(croots)[round(Im(croots), 10) == 0]
    ## the parametrization is for (x - x[i]), so need to shift the roots
    rroots <- rroots + x[i]
    ## real roots in (x[i], x[i + 1])
    xr[[i]] <- rroots[(rroots >= x[i]) & (rroots <= x[i + 1])]
    ## next piece
    i <- i + 1L
    }
  ## collapse list to atomic vector
  xr <- unlist(xr)
  ## make a plot?
  if (verbose) {
    curve(f, from = x[1], to = x[n_pieces + 1], xlab = "x", ylab = "f(x)")
    abline(h = y0, lty = 2)
    points(xr, rep.int(y0, length(xr)))
    }
  ## return roots
  xr
  }

It uses polyroot piecewise, first finding all roots on complex field, then retaining only real ones on the piecewise interval. This works because a cubic interpolation spline is just a number of piecewise cubic polynomials. My answer at How to save and load spline interpolation functions in R? has shown how to obtain piecewise polynomial coefficients, so using polyroot is straightforward.

Using the example data in the question, both RootSpline1 and RootSpline3 correctly identify all roots.

par(mfrow = c(1, 2))
RootSpline1(x, y, 2.85)
#[1] 3.495375 6.606465
RootSpline3(f3, 2.85)
#[1] 3.924512 6.435812 9.207171 9.886640