% 
clear;
a = -5; b = 5; 
n = 40; h = (b-a)/n;
xh = [a:h:b];
syms x 
fx = (x^2+1)^(-1);
dfx = diff(fx);
dfa = subs(dfx,'x',a);
dfb = subs(dfx,'x',b);
yh = subs(fx,'x',xh);
mu = ones(1,n-1)/2;
lambda = ones(1,n-1)/2;
[p,q] = chashang(xh,yh);
d = zeros(n+1,1);
d(2:n) = 6*p(1:n-1,4);
d(1) = 6*(p(1,3) - dfa)/h;
d(n+1) = 6*(dfb - p(n,3))/h;
A = diag(2*ones(1,n+1)) + diag([mu,1],-1) + diag([1,lambda],1);
M = A\d;
a3 = (M(2:n+1)-M(1:n))/(6*h);
a2 = M(1:n)/2;
a1 = (yh(2:n+1) - yh(1:n))'/h - h*(M(2:n+1) +2*M(1:n))/6;
a0 = yh(1:n)';


pp = spline(xh,[dfa yh dfb]);
[breaks,coefs,npolys,ncoefs,dim] = unmkpp(pp);
fprintf('\n err_a3=%f, err_a2=%f, err_a1=%f, err_a0=%f\n', ...
    norm(a3-coefs(:,1)), norm(a2-coefs(:,2)), ...
    norm(a1-coefs(:,3)), norm(a0-coefs(:,4)));

% plot the figures
np = 3;
n0 = np*n; h = (b-a)/n0;
xh0 = [a:h:b]; len = length(xh0);
xh =xh(1:n)';
for k = 1 : np
    index = [k:np:len];
    index = index(1:n);
    xht = xh0(index)'; 
    yht = a3.*(xht-xh).^3 + a2.*(xht-xh).^2 + a1.*(xht-xh) + a0;
    yh0(index) = yht;
end
yh0(n0+1) = yh(n+1);

yh1 = subs(fx,'x',xh0);
plot(xh0,yh1,'b.-',xh0,yh0,'ro-')