clear all
close all
nxy = 31;
xy = 2*(rand(2,nxy)-.5); vals = sum(xy.^2);
noisyvals = vals + (rand(size(vals))-.5)/5;
st = tpaps(xy,noisyvals);
xx = linspace(min(xy(1,:)),max(xy(1,:)),51);
yy = linspace(min(xy(2,:)),max(xy(2,:)),51);
[X,Y] = ndgrid(xx,yy);
mesh(X,Y,stval(st,{xx,yy}))
XY = [X(:)';Y(:)'];
DM = zeros(size(xy,2),size(XY,2));
for i = 1:2
  [xyd,XYd] = ndgrid(xy(i,:),XY(i,:));
  DM = DM+(XYd-xyd).^2;
end
% DM is the squared distance matrix
[centers,coefs] = stbrk(st); % get centers and coefs
n = size(centers,2);
vals = coefs(1:n)*phi(DM); % phi part
vals = vals+coefs(n+1)*XY(1,:)+coefs(n+2)*XY(2,:)+coefs(n+3); % linear part
figure
mesh(X,Y,reshape(vals,size(X)))