I want to solve the following differential equation
$$y'(t) = f(y(t),t)$$
where the right hand side is a function that has no closed form expression.
I have a vector of numerical values representing $f(y)$. It is time independent (each value corresponds to a grid point in y). How can I solve the differential equation?
I tried solving this as a loop as follows:
$$y(t+\Delta t) = y(t) + \Delta t f(y(t),t)$$
I use interpolation on the the second term. The problem is that this is very slow and I wonder if this can be solved more efficiently.
I would like to calculate how long it takes for $y(t_0)$ to reach a given value $y(T)$. For instance, in the graphs shown below, I want to know how long it takes to reach $\underline{y}$ for every point in the grid of $y$. How could I do this?

You mentioned Runge-Kutta methods in your tags to this question. Such methods assume the problem can be stated in the form $y'(t) = f(y(t),t).$ They do not assume that $f$ is sensitive to both of its inputs; you could have $f(y,t) = 1 + e^{-y}$ (or some other function that depended only on $y$) and that would be perfectly fine.
The standard application of a Runge-Kutta method to $y'(t) = f(y(t),t)$ means you choose a new value of $t$ at which to end each "step," so you end up with a list of function values $y(t_0)$ (where $t_0$ is your initial value of $t$), $y(t_1),$ $y(t_2),$ $y(t_3),$ and so forth. That is, you have a table of values of $y$ as a function of $t.$
Any of the usual Runge-Kutta methods should work and should produce such a table of values if you make sure to remember what the values of $y$ and $t$ were at the end of each step.
But it seems you really want $t$ as a function of $y.$ In fact that is what your "time to reach" figure represents if you consider it as the graph of a function with the input value on the horizontal axis and the output value on the vertical axis. If your $y$ is really just a real-valued strictly increasing function of time, as depicted in that figure, it is reasonable to think of the relation between $t$ and $y$ that way.
In that case, you could consider this equation: $$ \frac{dt}{dy} = \frac{1}{f(y)}. $$ Integrate both sides with respect to $y$: $$ t\big\rvert_{y=y_1} - t\big\rvert_{y=y_0} = \int_{y_0}^{y_1} \frac{1}{f(y)}\,dy. $$ If $y_0 = y(0),$ then $t\big\rvert_{y=y_0} = 0,$ so you can write $t(y)$ ($t$ as a function of $y$) as follows: $$ t(y) = \int_{y_0}^{y} \frac{1}{f(u)}\,du. $$ Your table of values of $f(y)$ lets you directly derive a table of values of $\frac{1}{f(u)},$ so now you're just integrating a function over an interval. You can use the trapezoid method, Simpson's method, or some other single-variable numeric integration, depending on what kind of guesses you want to make about the values of $\frac{1}{f(u)}$ in between the known values.
Runge-Kutta methods would still work, but seem to be overkill for a problem posed in that form.