function [a,b,c, res, Sres] = parabola(x,y)
% computes coefficients for best quadratic fit:
% y = a + b*x + c*x^2
% by solving the normal equations:
% 
m=length(x);
X1 = sum(x);
X2 = sum(x.*x);
X3 = sum(x.^3);
X4 = sum(x.^4);
Y1 = sum(y);
XY = sum(x.*y);
X2Y= sum(x.*x.*y);
AA = [m X1 X2; X1 X2 X3; X2 X3 X4];
bb = [Y1; XY; X2Y];
% Find polynomial coefficients for lest squares fit
xx= AA\bb;
a = xx(1);b = xx(2); c = xx(3);
% compute residual vector
res = a + b*x + c*x.*x - y;
Sres = sum(res.*res);
   