function [b, c, d] = cubic_spline_clamp(x, fun_val, df_val)

n       = length(x) - 1;
h       = zeros(n,  1);
A       = eye  (n+1,n+1);
rhs     = zeros(n+1,1);

for ii = 1:n
    h(ii) = x(ii+1) - x(ii);
end

A(1,1) = 2 * h(1);
A(1,2) = h(1);
rhs(1) = 3 * ( fun_val(2) - fun_val(1) ) / h(1) - 3 * df_val(1);
for ii  = 1:n-1
    A(ii+1,ii+1) = ( h(ii,1) + h(ii+1,1) ) * 2;
    A(ii+1,ii  ) = h(ii,1);
    A(ii+1,ii+2) = h(ii+1,1);
    rhs(ii+1,1)  = 3 * ( fun_val(ii+2,1) - fun_val(ii+1,1) ) / h(ii+1,1);
    rhs(ii+1,1)  = rhs(ii+1,1) - 3 * ( fun_val(ii+1,1) - fun_val(ii,1) ) / h(ii,1);
end
A(n+1,n  ) = h(n);
A(n+1,n+1) = 2 * h(n);
rhs(n+1)   = 3 * df_val(2) - 3 * ( fun_val(n+1) - fun_val(n) ) / h(n);

c = A \ rhs;

b = zeros(n,1);
d = zeros(n,1);
for ii = 1:n
    b(ii,1) = ( fun_val(ii+1,1) - fun_val(ii,1) ) / h(ii,1);
    b(ii,1) = b(ii,1) - h(ii,1) * ( c(ii+1,1) + 2 * c(ii,1) ) / 3;
    d(ii,1) = ( c(ii+1,1) - c(ii,1) ) / ( 3 * h(ii,1) );
end