function quadraticfit
   format compact
   m = 40;
   t = linspace(1,9,m);
   y = 1*ones(1,m) -2*t + (0.7)*t.^2  % start with an actual quadratic
   y = randompert(y,7)                % perturb the y-values a bit
   [alpha,beta,gamma] = lsq_linear_fit(t,y)

   f = @(t) alpha + beta*t + gamma*t.^2;
   hold off
   plot(t,y,'ro')
   range = max(t)-min(t);
   hold on
   myfplot(f,[min(t)-range/5,max(t)+range/5])
end

function [alpha,beta,gamma] = lsq_linear_fit(t,y) 
                         % Least squares straight line fit
                         % by solving the Normal Equations
   [q,m] = size(t);
   a = [ones(m,1) t' (t.^2)']
   b = y'
   [alpha,beta,gamma] = solve_normal_eqns(a,b)
end

function [alpha,beta,gamma]=solve_normal_eqns(a,b)
   xbar = ((a'*a)\(a'*b))';
   alpha = xbar(1);
   beta  = xbar(2);
   gamma = xbar(3);
end

function y=randompert(y,amp)
   y = y + amp*(rand(1,numel(y))-0.5);
end

function myfplot(f,interval)
	x = linspace(interval(1),interval(2),200);
	y = f(x);
	plot(x,y,'b');
end