n = 20;
x = [0, sort(rand(n-1,1))', 1];
f = @(x)real(1./(x-(.3+.1*1i))./(x-(.6-.2*1i)));
figure(1); clf;
plot(x,f(x),'rx','linewidth',2)
set(gca,'fontsize',16)

%%
h = x(2:n+1) - x(1:n);   % h_1 through h_N
inth = h(2:n-1);            % h_2 through h_N-1
sumh = h(1:n-1) + h(2:n);
A = diag(2*sumh) + diag(inth,-1) + diag(inth,1);
ff = f(x);

if (1),     % perturbation test
    ff(end-2) = 10;
    figure(1); clf; plot(x,ff,'rx','linewidth',2)
end

df = (ff(2:n+1) - ff(1:n))./h;       
b = 6*(df(2:n) - df(1:n-1))';
sigma = A \ b;

alpha = ff(2:n+1)./h - h.*[sigma' 0]/6;
beta = ff(1:n)./h - h.*[0 sigma']/6;

dx = 1/1e4;
xx = dx:dx:1-dx;
s = zeros(length(xx),1);
sigma = [0 sigma' 0];   % sigma(j+1) is old sigma(j)
for k = 1:length(xx)
    j = sum(xx(k)>x);
    hj = h(j);
    s(k) = (x(j+1)-xx(k))^3/(6*hj)*sigma(j) + ...
        (xx(k)-x(j))^3/(6*hj)*sigma(j+1) + ...
        alpha(j)*(xx(k)-x(j)) + beta(j)*(x(j+1)-xx(k));
end
figure(1); hold on;
plot(xx,s,'b-','linewidth',1);

%% 
plot(xx,f(xx),'r-','linewidth',1);