clear
format long
format compact

xl = -10.;
xr = 2.;
np = 200;

x = linspace(xl, xr, np);
y = airy(x);

plot(x, y, 'k-');
grid on
xlabel('x');
ylabel('Ai(x)');
title('Airy function and its zeroes');
hold on

help bisec
help secant
help newton

% By visual inspection conclude that the first (rightmost) negative root of
% Airy function is on the interval [-3., -1.], the second one is on the
% interval [-5., -3].

% First root

a1 = -3;
b1 = -1;

maxiter = 20;

x1b = bisec(@airy, a1, b1);
n1b = length(x1b)   % number of iterations/function evaluations
x1s = secant(@airy, a1, b1, maxiter);
n1s = length(x1s)   % number of iterations/function evaluations

airyprime = @(x) airy(1, x);  % use anonymous function of one variable for the derivative

x1n = newton(@airy, airyprime, a1, maxiter);
n1n = length(x1s)   % number of iterations
ninf = 2*n1n        % number of function evaluations

% are the roots consistent?
abs(x1b(end) - x1s(end))
abs(x1b(end) - x1n(end))
abs(x1s(end) - x1n(end))

% plot the root
plot(x1b(end), [0.], 'ro');

% Second root

a2 = -5;
b2 = -3;

x2b = bisec(@airy, a2, b2);
n2b = length(x2b)   % number of iterations/function evaluations
x2s = secant(@airy, a2, b2, maxiter);
n2s = length(x2s)   % number of iterations/function evaluations

x2n = newton(@airy, airyprime, a2, maxiter);
n2n = length(x2s)   % number of iterations
n2nf = 2*n2n        % number of function evaluations

% are the roots consistent?
abs(x2b(end) - x2s(end))
abs(x2b(end) - x2n(end))
abs(x2s(end) - x2n(end))
% plot the root
plot(x2n(end), [0.], 'ro');