function [c,R2,rout] = normalfit(x,y,basefun)
% Least-squares fit via solution to the normal equations
%
% Synopsis: 
% c = fitnorm(x,y,basefun)
% [c,R2] = fitnorm(x,y,basefun)
% [c,R2,r] = fitnorm(x,y,basefun)
%
% Input: x,y = vectors of data to be fit
% basefun = (string) m-file that computes matrix A with columns as
% values of basis basis functions evaluated at x data points.
%
% Output: c = vector of coefficients obtained from the fit
% R2 = (optional) adjusted coefficient of determination; 0 <= R2 <= 1
% r = (optional) residuals of the fit
if length(y)~= length(x)
    error('x and y are not compatible'); 
end

A = feval(basefun,x); % Coefficient matrix of overdetermined system
c = (A'*A)\(A'*y); % Solve normal equations, y(:) is always a column

if nargout>1
   r = y - A*c; % Residuals at data points used to obtain the fit
   [m,n] = size(A);
   R2 = 1 - (m-1)/(m-n-1)*(norm(r)/norm(y-mean(y)))^2;
   if nargout>2
     rout = r; 
   end
end