N = 20;
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);

plot(xfd, yfd, 'r.');
hold on
plot(xfd, yexact, '-');
grid on
