function MiddleRoot()
% Uses fzero to find the middle root of a cubic polynomial
% that has three real roots.
close all
a = [8 2 -5 1];
% Let p(x) = a(1) + a(2)x + a(3)x^2 + a(4)x^3
% Note that p'(x) = a(2) + 2a(3)x + 3a(4)x^2
% Let r1 and r2 be the zeros of p'(x). Because p has 3 real distinct
% roots, it follows that the middle root is in between r1 and r2
% and that p(r1)*p(r2)<0.

% Display p(x)...
xVals = linspace(-1.2,4.5);
plot(xVals,p(xVals,a),[-2 5],[0 0],'k')
shg
d = (2*a(3))^2 - 4*(3*a(4))*a(2);
% The roots of p'(x) are real and distinct
r1 = (-2*a(3) - sqrt(d))/(6*a(4));
r2 = (-2*a(3) + sqrt(d))/(6*a(4));
hold on
% Display the bracketing interval...
plot([r1 r2],[0 0],'r','linewidth',3)
plot([r1 r1],[-.4 .4],'r',[r2 r2],[-.4 .4],'r',r1,p(r1,a),'or',r2,p(r2,a),'or')

% The middle root is in between r1 and r2..
xStar = fzero(@(x) p(x,a),sort([r1 r2]));

plot(xStar,0,'.g','markersize',20)
title('p(x) = (x + 1)(x - 2)(x - 4) = 8 - 2x + 5x^2 + x^3','fontsize',16)
hold off

function z = p(x,a)
z = a(1) + a(2)*x + a(3)*x.^2 + a(4)*x.^3;