function [yp] = spline3clamped(xd,yd,xp,gammal,gammar)
% computes the values of S(x) at the points xp
% where S(x) is the clamped cubic spline interpolant
% for the data set (xd,yd) 
% satisfying the boundary conditions S'(x_1)=gammal, S'(x_n)=gammar

n=length(xd);

h=diff(xd); 
g=diff(yd)./h;

% set up n by n matrix M
M=diag([2*h(1); 2*(h(1:n-2)+h(2:n-1)); 2*h(n-1)], 0)...
    +diag(h,-1) + diag(h,1);

% set up rhs vector of length n
rhs=[3*(g(1)-gammal); 3*(g(2:n-1)-g(1:n-2)); 3*(gammar-g(n-1))];

% set up matrix of spline coefficients
S=zeros(n,4); % define the size
S(:,3)=M\rhs; % solve for c
for k=1:n-1
    S(k,1)=yd(k); % a
    S(k,4)=(S(k+1,3)-S(k,3))/(3*h(k)); % d
    S(k,2)=g(k) - (2*S(k,3)+S(k+1,3))*h(k)/3; % b
end

Srev=S(1:n-1,4:-1:1); % discard the last row, reverse order for mkpp
pp=mkpp(xd,Srev);
yp=ppval(pp,xp);

end

