% interpolated points
P = [0 1; 1 3; 6 2; 3 0];
plot([P(:,1);P(1,1)],[P(:,2);P(1,2)],'*b-')
hold on

M = toeplitz([1/2 0 0 1/2],[1/2 1/2 0 0]);
L = toeplitz([3/4 1/8 0 1/8]);

% solution to the interpolation problem
cP = pinv(M)*P;
tP = L\P;

% subdivision curve for each solution
stP = SubdivS(tP,5,[1/4 3/4 3/4 1/4]);
scP = SubdivS(cP,5,[1/4 3/4 3/4 1/4]);

plot([stP(:,1);stP(1,1)],[stP(:,2);stP(1,2)],'*g-')
plot([scP(:,1);scP(1,1)],[scP(:,2);scP(1,2)],'*m-')