function prog1
    % This function is really just a script - it calls the secant method
    % code (found below) and generates all of the tables and plots
    % included in my solution.
    %
    % Note: the generated tables are in LaTeX format.

    % exercise 1

    f = 'sin(x)/x^2 - cos(x)/x';

    % compute roots
    
    x1 = secant(f, 3, 3.1, eps, eps, 100 );
    x2 = secant(f, 7, 7.1, eps, eps, 100 );
    x3 = secant(f, .1, .9, eps, eps, 100 );

    % generate table
    
    T = [];
    T(2,1:length(x1)) = x1;
    T(3,1:length(x2)) = x2;
    T(4,1:length(x3)) = x3;
    T(1,:) = (1:size(T,2)) - 1;
    fprintf('$p_%d$ & %.4f & %.4f & %.4e \\\\ \n', T);
    
    % exercise 2
    
    f = '(x - 1)^4 * (x + 1)';
    
    % compute roots
    
    x1 = secant(f, .5, .6, eps, eps, 100 );
    x2 = secant(f, -1.5, -1.4, eps, eps, 100 );
    
    % generate tables
    
    e1 = abs(x1 - 1);
    r1 = e1(2:end) ./ e1(1:end-1);
    fprintf('%d & %.4f & %.4e & %.4e \\\\\n', [0:length(x1)-1; x1; e1; [r1 0]]);
    
    e2 = abs(x2 + 1);
    r2 = e2(2:end) ./ e2(1:end-1) .^ ((1 + sqrt(5))/2);
    fprintf('%d & %.4f & %.4e & %.4e \\\\\n', [0:length(x2)-1; x2; e2; [r2 0]]);
    
    % exercise 3
    
    p = '(x - 2)^5';
    q = '(x - 2)^5 + 1e-6 * x^3';
    
    % make perturbation plot
    
    figure(1);
    t = linspace(1.85, 2.15, 1000);
    P = inline(vectorize(p));
    Q = inline(vectorize(q));
    plot(t, 0*t, 'g-', t, P(t), 'b-', t, Q(t), 'r-');
    axis('tight');
    legend('x axis', '(x - 2)^5', '(x - 2)^5 + 10^{-6} x^3', 'Location', 'SouthEast');
    
    % compute new root
    
    x = secant(q, 2.1, 2, eps, eps, 100);
    r = x(end)
    ratio = (r - 2) / 1e-6
    
    % exercise 4
    
    p = '(x - 18)*(x-19)*(x-20)*(x-21)*(x-22)';
    q = '(x - 18)*(x-19)*(x-20)*(x-21)*(x-22) + 6e-7 * x^5';

    % make perturbation plot
    
    figure(2)
    P = inline(vectorize(p));
    Q = inline(vectorize(q));
    t = linspace(17.8, 22.2, 1000);
    plot(t, 0*t, 'g-', t, P(t), 'b-', t, Q(t), 'r-');
    axis('tight');
    legend('x axis', 'p(x)', 'p(x) + 6 \times 10^{-7} x^5');
    
    % compute new roots
    
    x1 = secant(q, 18, 18.1, eps, eps, 100);
    r1 = x1(end)
    e1 = abs(r1 - 18)/18
    
    x2 = secant(q, 21, 21.1, eps, eps, 100);
    r2 = x2(end)
    e2 = abs(r2 - 21)/21
    
    x3 = secant(q, 22, 21.9, eps, eps, 100);
    r3 = x3(end)
    e3 = abs(r2 - 22)/22
end

function [x, y, err, i] = secant(f, p0, p1, xtol, ytol, maxit) 
    F = inline(f);

    x(1) = p0;
    y(1) = F(p0);
    
    x(2) = p1;
    y(2) = F(p1);
    
    err = x(2) - x(1);
    
    i = 2;
  
    while i < maxit && abs(err) > xtol && abs(y(i)) > ytol
        i = i + 1;
        x(i) = x(i-1) - y(i-1) * (x(i-1) - x(i-2)) / (y(i-1) - y(i-2));
        y(i) = F(x(i));
        err = x(i) - x(i-1);
    end
end