function y=tps(r)% This is the thin-plate splineif r < 0.000000000001 y=0;else y=r^2*log(r);endfunction y=fun(point)% my target functionx=point(1);z=point(2);y=7-4*x^2+z^3;function y=interpmat(points)% this file computes the interpolation matrix[n,m]=size(points);for i=1:nipoint=points(i,:);for j=1:njpoint=points(j,:);mat(i,j)=tps(norm(ipoint-jpoint));endmat(i,n+1)=1;mat(n+1,i)=1;mat(i,n+2)=ipoint(1);mat(n+2,i)=ipoint(1);mat(i,n+3)=ipoint(2);mat(n+3,i)=ipoint(2);endy=mat;function y=righthandside(points)% find the right hand side of the interpolation system[n,m]=size(points);for i=1:ny(i)=fun(points(i,:));endy(n+1)=0;y(n+2)=0;y(n+3)=0;function y=testinterp(n,m)% this file test our interpolation routinepoints=rand(n,2);% compute interpolation matrixmat=interpmat(points);% create the right hand siderhs=righthandside(points);solution=inv(mat)*rhs';% test solutiontpoints=rand(m,2)maxerr=testsol(points,tpoints,solution);