%% Find Critical Points and Normal Lines
 
% Critical points
CP = sort(roots(p), 'ascend');
plot(CP(1)+50, 0, 'bo', 'LineWidth', 2, 'MarkerSize', 10); % top critical point
hold on
plot(CP(2)+50, 0, 'mo', 'LineWidth', 2, 'MarkerSize', 10); % middle critical point
hold on
plot(CP(3)+50, 0, 'ro', 'LineWidth', 2, 'MarkerSize', 10); % bottom critical point
 
% Slope of polynomial at critical points
x2 = linspace(CP(1)+45, CP(1)+55, 2);
f3 = polyval(p, x2);
slope1 = 10/(f3(2)-f3(1));
x3 = linspace(CP(2)+45, CP(2)+55, 2);
f4 = polyval(p, x3);
slope2 = 10/(f4(2)-f4(1));
x4 = linspace(CP(3)+45, CP(3)+55, 2);
f5 = polyval(p, x4);
slope3 = 10/(f5(2)-f5(1));
 
% Top intersection point and normals
m1 = -(slope1)^-1;
m2 = -(slope2)^-1;
m3 = -(slope3)^-1;
b1 = CP(1)+50;
b2 = CP(2)+50;
b3 = CP(3)+50;
 
A = [m2, 1; m1, 1];
b = [b2; b1];
topIP = A\b; % intersection point of normals
plot(topIP(2), -topIP(1), 'bo', 'LineWidth', 2, 'MarkerSize', 5); 
hold on
 
if topIP(1) < 0
    x5 = linspace(0, -topIP(1), -topIP(1));
else
    x5 = linspace(-topIP(1), 0, topIP(1));
end
norm1 = m1*x5 + b1;
plot(norm1, x5, 'b'); % bottom line
hold on
norm2 = m2*x5 + b2;
plot(norm2, x5, 'b'); % top line
hold on
 
% Bottom intersection point and normals
A = [m3, 1; m2, 1];
b = [b3; b2];
botIP = A\b; % intersection point of normals
plot(botIP(2), -botIP(1), 'ro', 'LineWidth', 2, 'MarkerSize', 5); 
hold on
 
if botIP(1) < 0
    x6 = linspace(0, -botIP(1), -botIP(1));
else
    x6 = linspace(-botIP(1), 0, botIP(1));
end
norm2 = m2*x6 + b2;
plot(norm2, x6, 'r'); % bottom line
hold on
norm3 = m3*x6 + b3;
plot(norm3, x6, 'r'); % top line
hold on
 
camroll(-90);
