function [a b c]=lsquare(x,y,n)
  for i=1:n+1
    for j=1:n+1
      a(i,j)=sum(x.^(j+i-2));
    endfor
  endfor
  for i=1:n+1
    b(i)=sum(y.*x.^(i-1));
  endfor
  b=b';
  c=a\b;
  lsquarev=[a b c];
endfunction