Here is what I meant:

```
function out=f(par,xij,yij,x,y)
i = par(1); j = par(2);
out = [xij-interp2d(i,j,1:10,1:9,sx,'natural')
yij-interp2d(i,j,1:10,1:9,sy,'natural')];
endfunction
sx = splin2d(1:10,1:9, x,'natural');
sy = splin2d(1:10,1:9, y,'natural');
sP = splin2d(1:10,1:9, P,'natural');
// goal: get a “P-Value” for point k
xk = 0;
yk = 9e6;
// determine parameters value for given (xk,yk)
par = fsolve([5;5],list(f,xk,yk,x,y))
// compute P-value for computed parameters (i,j)
pk = interp2d(par(1), par(2), 1:10, 1:9, sP, 'natural')
```

Compared to the bilinear interpolation, the 2d spline catches more curvature information and should (or should not) better extrapolate. The default type type of spline (not a knot) does not give good extrapolation that’s why ‘natural’ is used (it has a different meaning than the ‘natural’ of interp2d).

But depending of your goal, maybe interpolation is not the good way to go. I you want to extrapolate, you should rather use a global model (e.g. polynomial) with a small number of parameters.

Here are two illustrations of the extrapolation (mesh is in red) depending of the used interpolation type:

```
// bilinear interpolation
[i,j]=meshgrid(0:0.5:11,0:0.5:10);
X=linear_interpn(i,j,1:10,1:9,x,'natural');
Y=linear_interpn(i,j,1:10,1:9,y,'natural');
P=linear_interpn(i,j,1:10,1:9,P,'natural');
k = 3:$-2;
mesh(X,Y,P,'edgecolor','red')
mesh(X(k,k),Y(k,k),P(k,k))
```

```
// 2d spline interpolation
[i,j]=meshgrid(0:0.5:11,0:0.5:10);
X=interp2d(i,j,1:10,1:9,sx,'natural');
Y=interp2d(i,j,1:10,1:9,sy,'natural');
P=interp2d(i,j,1:10,1:9,sP,'natural');
k = 3:$-2;
mesh(X,Y,P,'edgecolor','red')
mesh(X(k,k),Y(k,k),P(k,k))
```

I don’t like the amplification of oscillations in the lower right corner…

S.