function findmaxima

data=[
    1  3.0  1.0 
    2  5.0  0.8
    3  7.0  0.8
    4  6.5  0.7
    5  3.0  1.0
];

x = data(:,1);  % 1st column
y = data(:,2);  % 2nd column
s = data(:,3);  % 3rd column


%
% number of data points
%
n=length(x);


%
% Form the M matrix
%
for r=1:3
    for c=1:3
        M(r,c) = sum((x.^(r+c-2))./(s.^2));
    end
end


%
% Form the N vector
%
for r=1:3
    N(r,1) = sum((x.^(r-1)).*(y)./(s.^2));
end


%
% Print out the results. C is the inverse of M, and forms variances and
% covariances.
%
M
N
C = inv(M)
A = C*N
a = A(1);
b = A(2);
c = A(3);



%
% Find the Chi-square of the fit
%
chi2 = sum((y-(a + b*x + c*x.^2)).^2 ./ s.^2)
chi2_dof = chi2 / (n-3)


%
% Prepare the fitted polynomial to be drawn
%
dx = (max(x)-min(x))/1000;
xx=min(x)-dx*100:dx:max(x)+dx*100;
yy = a + b*xx + c*xx.^2;
ss = sqrt(C(1,1) + C(1,2)*xx + (C(2,2)+2*C(1,3))*xx.^2 + (2*C(2,3))*xx.^3 + (C(3,3))*xx.^4);

%
% Plot the data and fitted-polynomial on top of it
%
hold off
errorbar(x,y,s,'ro')
hold on
plot(xx,yy)
plot(xx,yy+ss,'k:')
plot(xx,yy-ss,'k:')
hold off


%
% Find the maxima point and its error
%
xmax = -b/(2*c)
xmax_err = sqrt((1/(2*c))^2*C(2,2) + (b/(2*c^2))^2*C(3,3) - 2*(b/(4*c^3))*C(2,3))


%
% Mark the maxima region
%
hold on
ax = axis;
plot([xmax xmax],[ax(3) ax(4)],'--k')
plot([xmax-xmax_err xmax-xmax_err],[ax(3) ax(4)],':k')
plot([xmax+xmax_err xmax+xmax_err],[ax(3) ax(4)],':k')
hold off
h = text(xmax+1.5*xmax_err,ax(3)+0.1*diff(ax(3:4)),sprintf('x_{max} = %.3f \\pm %.3f',xmax, xmax_err));
set(h,'FontSize',16)



%
% Creates a pdf file of the plot
%
print -dpdf findmaxima.pdf
