format long;
% The function
f = @(x) (6425*x.^8 -12012*x.^6 + 6930*x.^4 - 1260*x.^2 + 35)/128;
% its derivative
fp = @(x) (8*6425*x.^7 -6*12012*x.^5 + 4*6930*x.^3 - 2*1260*x + 35)/128;

x = -1:0.01:1;
plot(x, x*0, 'k-');
hold on;
plot(x,f(x));
hold on;

N = 3;  % Number of bisection steps

% The numbers in the interval array are chosen by eyeballing the minimums
% and the maximums, after plotting the data in the interval of interest.
intervals = [-1 -0.9 -0.66 -0.36 0 0.36 0.66 0.9 1];

answers = zeros(size(intervals,2)-1,1); % An array for keeping the zeros.
steps = zeros(size(intervals,2)-1,2);   % A matrix for keeping the number of steps.
                                        % First column for bisection steps
                                        % and secton column for Newton
                                        % steps.
tol = 1e-15;    % Turns out eps is a bit too small.
                % So I changed the tolerance to something slightly larger.
                % This is used for terminating Newton's steps.
relErr = zeros(size(intervals,2)-1,1);

for i=1:size(intervals,2)-1
    a = intervals(i);
    b = intervals(i+1);
    j = 1;
    % N bisection steps
    while j < N
        newPoint = (a+b)/2;
        if f(a)*f(newPoint) > 0
            oldPoint = a;
            a = newPoint;
        else
            oldPoint = b;
            b = newPoint;
        end
        j = j + 1;
    end
    
    k = 1;
    % Newton steps until we have to terminate.
    while abs(oldPoint-newPoint) > tol*abs(oldPoint)
        oldPoint = newPoint;
        newPoint = newPoint - f(newPoint)/fp(newPoint);
        k = k + 1;
        % disp(abs(oldPoint-newPoint));
    end
    relErr(i) = abs((oldPoint-newPoint)/oldPoint);
    answers(i) = newPoint;
    steps(i,1) = j;
    steps(i,2) = k;
end
disp('  Answer                Relative Error      Bisection           Newton          Total Steps')
disp([answers relErr steps(:,1) steps(:,2)  steps(:,1)+steps(:,2)]);
plot(answers, f(answers), 'or');