interpolating scattered points using radial basis function

We suppose that we have sampled the value of a function

Now i want to include a linear term in the interpolating function g(x)

.

But with 3 new unknows i need 3 new equations which are as follow:

I have modified the linear system (1) to include the new unkowns and equations

so i have tried to solve the modifed problem with the surface data provided in samples.mat but end up getting an errors saying 'Arrays have incompatible sizes for this operation.'. I dont know whats causing the error or how to fix it. Any assistance would be appreciated

load samples.mat %samples.mat outputs pts, f, c_provided(this can be ignored)
phi = @(r) r.^2 .* log(r);  % thin plate spline RBF
x1=pts(1:1089,1);
x2=pts(1:1089,2);
A=phi(pdist2(pts, pts));
for i=1:length(x1);
    A(i,i)=0; %fixing the diagonals that equal NaN
end
%iteration to modify the matrix
for i=1:1089;
    A(i,1090)=1;
    A(i,1091)=x1(i);
    A(i,1092)=x2(i);
    A(1090,i)=1;
    A(1091,i)=x1(i);
    A(1092,i)=x2(i);
end

f(1090:1092,1)=0;
c=A\f;
a_0=c(1090,1);
a_1=c(1091,1);
a_2=c(1092,1);
g = @(x) phi(pdist2(x, pts)).*c_provided + a_0 + a_1.*x1 +a_2.*x2; %mathlab say the error is here
[X1, X2] = meshgrid(linspace(0,1,gridsize));
G = g([X1(:) X2(:)]);
G = reshape(G, size(X1));
figure
surf(X1, X2, G)
axis vis3d
xlabel x_1, ylabel x_2
zlabel f