function [A,B,x,y,z]=mquad(x_cen,y_cen)
% Produces Hardy multiquadric interpolant to the peaks function,
% given centers with xy coordinates in x_cen and y_cen.
%
ydata=peaks(x_cen,y_cen);
ydata=[ydata;0];
x=x_cen;
y=y_cen;
N=length(x_cen);
[XX,YY]=rbf_inputs(x_cen,y_cen);
A=(1+XX.^2+YY.^2).^(1/2);
B=ones(N+1,N+1);
B(1:N,1:N)=A;
B(N+1,N+1)=0;
coefs=pinv(B)*ydata;
x_val=-3:0.1:3;
y_val=-3:0.1:3;
M=length(x_val);
[x,y]=meshgrid(x_val,y_val);
[xx,yy]=rbf_inputs(x_cen,y_cen,x(:),y(:));
z=reshape(coefs(1:N)'*(1+xx.^2+yy.^2).^(1/2) + coefs(N+1),M,M);
subplot(2,1,1), peaks(x,y)
subplot(2,1,2), surf(x,y,z), axis([-3 3 -3 3 -8 8]),...
    title('Multiquadric Interpolant')



