This is very helpful. Definitely binary-search would be much better for large N. In my case, N was about 10 with most interpolations being done in the first 2-3 and so binary-search didn't really help. But, in general, of course it's much better.
In my use-case, also scipy.interpolate.interp1d doesn't do what I needed for the situation where yvals.shape = (len(xnew), len(xvals)) with both xvals and xnew being 1-d. I needed the output to be 1-d with length the same as xnew. But, the output of scipy.interpolate.interp1d would have shape (yvals.shape, len(xnew)) --- i.e. an cross-product instead of a "zip" operation.