Allow cubic interpolant for 2D data

[TomTranter]

Description

Currently the scipy interpolate.interp2d function is called if there are 2 x arguments and y is 2d but we only allow "linear" interpolants. The function will accept "cubic" and "quintic" also. It would be useful to have these options accessible. Were they prevented for any good reason? Convergence issues maybe?

Motivation

Functions that interpolate 2 variables will be more accurately represented.

Possible Implementation

No response

Additional context

No response

[valentinsulzer]

Don't think they were prevented for any particular reason, there just wasn't a good reason to have them.
We should definitely add them. But we should be wary of the fact that they don't preserve monotonicity

[rtimms]

Fixed in #2229