mq2d(x,xc,c)

function [A,Ax,Ay,Axx,Ayy] = mq2d(x,xc,c)
% 2-D multiquadric radial basis function
% f = @(r,c) sqrt((c*r).^2 + 1);

f = @(r,c) sqrt((c*r).^2 + 1);
N = size(x,1);

dx = x - repmat(xc,N,1);
r = sqrt(sum(dx.^2,2));

% Interpolation
A = f(r,c);

if nargout > 1
% 1-st derivative w.r.t x
Ax = c^2*dx(:,1)./A;
    if nargout > 2
    % 1-st derivative w.r.t y
    Ay = c^2*dx(:,2)./A;
        if nargout > 3
        % 2-nd derivative w.r.t x
        Axx = c^2*(1 + (c*dx(:,2)).^2).*(A.^-3);
            if nargout > 4
            % 2-nd derivative w.r.t y
            Ayy = c^2*(1 + (c*dx(:,1)).^2).*(A.^-3);
            end
        end
    end
end