3D Lagrange Interpolation

Hello, I'm trying to interpolate a polynomial function through a set of points obtained from a simulation. However I'm always getting the error "Z must be a matrix, not a scalar or vector" even though Z is actually a Matrix. And also getting the error:

in lagrange2D (line 43)

surf(x(1):(x(end)-x(1))/(size(A,1)-1):x(end),y(1):(y(end)-y(1))/(size(A,1)-1):y(end),A);

I have this data

%just xh for example, but there is also yh and F1h, and F1h = f(xh,yh)
xh =

  Columns 1 through 9

    0.0112    0.0112    0.0112    0.0112    0.0112    0.0112    0.0112    0.0090    0.0090

  Columns 10 through 18

    0.0090    0.0090    0.0090    0.0090    0.0090    0.0140    0.0140    0.0140    0.0140

  Columns 19 through 27

    0.0140    0.0140    0.0140    0.0135    0.0135    0.0135    0.0135    0.0135    0.0135

  Columns 28 through 36

    0.0135    0.0075    0.0075    0.0075    0.0075    0.0075    0.0075    0.0075    0.0064

  Columns 37 through 45

    0.0064    0.0064    0.0064    0.0064    0.0064    0.0064    0.0106    0.0106    0.0106

  Columns 46 through 49

    0.0106    0.0106    0.0106    0.0106

and then the function itself

    function lagrange2D
%% STARTVARIABLES
N = 15;                                           
xa = linspace(min(xh), max(xh), N);                   
ya = linspace(min(yh), max(yh), N);                  
[X,Y] = meshgrid(xa, ya);                           
Z = griddata(xh, yh, F1h, X, Y);
step = 0.2;
%% INTERPOLATION
% erst kurven in x richtung ausrechnen
xCurves={};
for i=1:size(X,1)
    % x vector must be a column vector for the lagrange function
    x = X(i,:)';
    % z vector must be a column vector for the lagrange function
    z = Z(i,:)';
    p=[];
    for j = x(1):step:x(end)
        % interpolate for every parameter value j and add it to p
        p = [p,lagrange(x,z,j)];
    end
    % save curves in x direction
    xCurves{i} = p;
end
y = Y(:,1);
% matrix for the graphical outpu
A=[];
% interpolate in y-direction
for i=1:length(xCurves{1})
    p=[];
    z=[];
    for l=1:length(y)
        z = [z;xCurves{l}(i)];
    end
    for j = y(1):step:y(end)
         % interpolate for every parameter value j and add it to p
        p = [p;lagrange(y,z,j)];
    end
    A = [A,p];
end
%% GRAPHICAL OUTPUT
% plot surface
surf(x(1):(x(end)-x(1))/(size(A,1)-1):x(end),y(1):(y(end)-y(1))/(size(A,1)-1):y(end),A);
hold on;
% plot points
%for i=1:size(X,1)
    %for j=1:size(Y,1)
        %p = plot3(X(i,j),Y(i,j),Z(i,j));
        %set(p,'Marker','.');
        %set(p,'MarkerSize',30);
    %end
%end
    end

Any thoughts on what am I doing wrong in here?

Thanks in Advance

[Eduardo Hulse]

I forgot to mention that my lagrange(x,z,t) used has the script

%% nested functions
function res = lagrange(x,z,t)
% lagrange interpolation in 1d
% input: x...vektor, welcher die x (bzw. y) Werte beinhaltet
res = 0;
for i=1:length(x)
    % delete the component not needed
    temp = [x(1:i-1);x(i+1:end)];
    
    denominator = x(i)-temp;
    
    numerator = t-temp;
    res = res + z(i)*(prod(numerator))/(prod(denominator));
end
end

[John D'Errico]

Could I point out that Lagrange interpolation is a nasty, ill-tempered thing to do? Even Lagrange would not be using that method if he had a choice. And the fact is, you don't need to use it, as there are far better methods available. But to make things worse, you want to use that in three dimensions? Are you serious? This is just a flat out bad idea.

[Eduardo Hulse]

I know it is a bad idea and that is one of the purposes. I will test the same data to interpolate a function using Kriging method and then compare it with lagrange.