function [p] = my_ls(x,y,m,n);

p=zeros(n+1,1);  

% Form the coefficient matrix. Input x and y should be colomn vecors of m+1

A=zeros(m+1,n+1);

for i=1:m+1,
  for j=1:n+1
     A(i,j) = x(i)^(j-1); 
  end
end

% Form the normal equations

 B = A'*A;  b=A'*y; 

 p = B\b;

end