function b = trend(X,m,lab,ngrid)
if nargin < 4, ngrid = 10; end
if nargin < 3, lab = 0; end
[N,nc] = size(X);  % determine no of data (rows)
X1 = polygen(X(:,1:2), m);  % generate polynomial terms
X1 = [ones(N,1) X1 X(:,3)]; % expand by adding a first column of ones
S = X1'*X1;          % calculate sums of squares and products matrix
[n,n2] = size(S);
SS = S(1:n-1,1:n-1);   % normal equations are  SS*b = y 
y = S(1:n-1,n);
b = SS\y;              % solution of normal equations
yhat = X1(:,1:n-1)*b;  % col. vector of expected values (y hat)
e = X1(:,n) - yhat;    % col. vector of residuals
T = cotomat(t1,t2);  % generate list of coords
Xt = polygen(T, m);  % generate polynomial terms
[rt ct] = size(Xt);
Xt = [ones(rt,1) Xt];
that = Xt*b;         % trend at each grid point
Y = dattogr([T that], ngrid+1, ngrid+1);
v = contvec(yhmin,yhmax);  
c = contour(Y,v);
if lab == 1
   clabel(c);
end
axis image;  % change to axis('equal') for MATLAB 4
s = num2str(m);
title(['Contours on order ',s,' Trend']);
xlabel('x'), ylabel('y');
