nmin = 4;
nmax = 12;
np = nmax - nmin + 1;

err = zeros(np,1);
sz = zeros(np,1);

for i = nmin:nmax
    N = 2^i;
    sz(i-nmin+1) = N;
    err(i-nmin+1) = myfderr(N);
end

loglog(sz, err, 'r.-');
grid on
hold on
loglog(sz, sz.^(-2)/2^nmax, '--');
legend('Numerical experiment', 'N^{-2} guide line')

function err = myfderr(N)
    a = 1;
    b = 2;
    alpha = 1;
    beta = 2;

    xfd = linspace(a, b, N+2);
    delta = xfd(2) - xfd(1);

myp = @(x) -2 ./ x;
myq = @(x)  2 ./ x.^2;
myr = @(x)  sin(log(x)) ./ x.^2;

    p = myp(xfd(2:(N+1)));
    q = myq(xfd(2:(N+1)));
    r = myr(xfd(2:(N+1)));

    f = delta^2 * r;
    f(1) = f(1) - alpha*(1 + p(1)*delta/2);
    f(N) = f(N) - beta*(1 - p(N)*delta/2);

    d1 = 1 + delta*p(2:N)/2;
    d1 = [d1 0];
    d2 = -(2 + delta^2 * q);
    d3 = 1 - delta*p(1:N-1)/2;
    d3 = [0 d3];
    A = spdiags([d1' d2' d3'], [-1 0 1], N, N);

    yfd = A \ f';
    yfd = [alpha yfd' beta];
    yexact = myexact(xfd);

    err = norm(yfd - yexact, inf);
end

