clc, clear, close all;


function [root,errors] = bisectionmethod(xl,xu,tol,m_iters,f)
if (f(xl)*f(xu) >= 0)
    root = NaN;
    errors = [];
    return
end

errors = zeros(1,m_iters);

cp = (xl+xu)/2;

for i = 1:m_iters
    c = (xl+xu)/2;
    errors(i) = abs((c-cp)/cp);
    if abs(f(c)) < tol || abs(xu-xl)/2 < tol
        root = c;
        fprintf("Converged in %d iterations\n",i);
        errors = errors(1:i);
        return;
    end
    if (f(xl) * f(c) < 0)
        xu = c;
    else
        xl = c;
    end
end

if i == m_iters
    disp("Max Iterations reached.")
    root = c;
end

end

function [root,errors] = newtonraphsonmethod(x0,tol,m_iters,f,df)
c = x0;

errors = zeros(1,m_iters);

for i = 1:m_iters
    if (abs(df(c)) < 1e-12)
        error('Division By Zero.')
    end
    cn = c - f(c)/df(c);
    errors(i) = abs((cn-c)/cn);
    if abs(cn-c)< tol
        c = cn;
        root = c;
        fprintf("Converged in %d iterations\n",i);
        errors = errors(1:i);
        return;
    end
    c = cn;
end

if i == m_iters
   disp("Max Iterations reached.")
   root = c;
end

end

function [root,errors] = secantmethod(x0,x1,tol,m_iters,f)
c = x0;
cn = x1;

errors = zeros(1,m_iters);

for i = 1:m_iters
    if abs(f(cn)-f(c)) < 1e-12
        error('Division By Zero.')
    end
    cnn = cn - f(cn) * (cn-c)/(f(cn)-f(c));
    errors(i) = abs((cnn-cn)/cnn);
    if abs(cnn - cn) < tol
        root = cnn;
        fprintf("Converged in %d iterations\n",i);
        errors = errors(1:i);
        return;
    end
    c = cn;
    cn = cnn;
end

if i == m_iters
   disp("Max Iterations reached.")
   root = cnn;
end
end

syms x f(x) d(x)

%defined from (0,+inf)
f(x) = 4 * x^9 + 8 * x^3 - log(5 * x);

d(x) = diff(f(x));

f = matlabFunction(f(x));
df = matlabFunction(d(x));

tol = 1e-6;
m_iters = 50;

% bisection
[root1bs, errors1bs] = bisectionmethod(0.2, 0.25, tol, m_iters, f);
[root2bs, errors2bs] = bisectionmethod(0.45, 0.5, tol, m_iters, f);

% nr method
[root1nr, errors1nr] = newtonraphsonmethod(0.2, tol, m_iters, f, df);
[root2nr, errors2nr] = newtonraphsonmethod(0.5, tol, m_iters, f, df);

% secant method
[root1sc, errors1sc] = secantmethod(0.15, 0.2, tol, m_iters, f);
[root2sc, errors2sc] = secantmethod(0.4, 0.45, tol, m_iters, f);

fstring = "Root 1 BS %f, Root 2 BS %f\nRoot 1 NR %f, Root 2 NR %f\nRoot 1 SC %f, Root 2 SC %f";
fprintf(fstring,root1bs,root2bs,root1nr,root2nr,root1sc,root2sc)

figure;
plot(1:length(errors1bs), errors1bs, '-o','LineWidth',2);
title('Bisection Method - Approximate Error for Root 1 (0.2, 0.25)');
xlabel('Iterations');
ylabel('Approximate Relative Error');
saveas(gcf, 'Bisection_Root1_Error.png');
grid on;

figure;
plot(1:length(errors2bs), errors2bs, '-o','LineWidth',2);
title('Bisection Method - Approximate Error for Root 2 (0.45, 0.5)');
xlabel('Iterations');
ylabel('Approximate Relative Error');
grid on;
saveas(gcf, 'Bisection_Root2_Error.png');

figure;
plot(1:length(errors1nr), errors1nr, '-o','LineWidth',2);
title('Newton-Raphson Method - Approximate Error for Root 1 (x0 = 0.2)');
xlabel('Iterations');
ylabel('Approximate Relative Error');
grid on;
saveas(gcf, 'Newton_Raphson1_Error.png');

figure;
plot(1:length(errors2nr), errors2nr, '-o','LineWidth',2);
title('Newton-Raphson Method - Approximate Error for Root 2 (x0 = 0.5)');
xlabel('Iterations');
ylabel('Approximate Relative Error');
grid on;
saveas(gcf, 'Newton_Raphson2_Error.png');

figure;
plot(1:length(errors1sc), errors1sc, '-o','LineWidth',2);
title('Secant Method - Approximate Error for Root 1 (x0 = 0.15, x1 = 0.2)');
xlabel('Iterations');
ylabel('Approximate Relative Error');
grid on;
saveas(gcf, 'Secant1_Error.png');

figure;
plot(1:length(errors2sc), errors2sc, '-o','LineWidth',2);
title('Secant Method - Approximate Error for Root 2 (x0 = 0.4, x1 = 0.45)');
xlabel('Iterations');
ylabel('Approximate Relative Error');
grid on;
saveas(gcf, 'Secant2_Error.png');
