clear all
close all
%%
%%Defining functions
%P8
P=(@(x) (6435*x.^8-12012*x.^6+6930*x.^4-1260*x.^2+35)/128);

%The derivitive of P8
dP=(@(x) (51480*x.^7-72072*x.^5+27720*x.^3-2520*x)/128);

%%
n=10; % of regions which the x axis is divided into initialy
%Must have only one root of P8 in each region.
%Inspection of the graph shows that n=10 is enough.
regionends=[-1:2/n:1];%the endpoints of the regions

%%

% test plot
X=-1:0.01:1;
figure
plot(X,P(X),'LineWidth',2)
hold on
plot([1 -1],[0 0],'k')
set(gca,'XTick',regionends)% x ticks mark out the regions
%plot(X,0.1*dP(X),'r')
xlabel('$x$','interpreter','latex','FontSize', 18)
ylabel('$P_8(x)$','interpreter','latex', 'FontSize', 18)

%%
%Root solving
maxiter=1000;%Max iterations for Newton's method
tolerance=1e-12;%tolarence for Newton's method
roots=nan*[1:1:n];%Blank list to load roots into 
%(will use isnan to remove unfilled elements)
for i=1:n % loop over regions
    if P(regionends(i))*P(regionends(i+1))<0% if P has a difrent sign 
        %on each end of a region then that region contains a root.
        
        %Bisection algorithm
        XL=regionends(i);
        XR=regionends(i+1);
        for ii=1:3%three iterations of the Bisection algorithm
            xmiddle=(XL+XR)/2;
            if P(XL)*P(xmiddle)<0; XR=xmiddle;else XL=xmiddle;end;
        end
        
        %Newton's Method
        x=(XL+XR)/2;% initial guess for Newton's method from
        %bisection method results.
        iter=1;%reset iteration counter
        accuracy=1;%reset accuracy
        while accuracy>tolerance && iter<maxiter;
            x=x-feval(P,x)/feval(dP,x);
            accuracy=abs(feval(P,x));
        end
       
        roots(i)=x;   %load root into list   
        
    end
      
end

roots=roots(isnan(roots)==0)% removes NaNs from the vector roots
% and print roots to command window.

%plot roots on graph 
plot(roots,0*roots,'ok','LineWidth',2,'MarkerFaceColor','r')
print( '-depsc2', 'Q1fig' )%save figure for report