function Z = Zero(P,T)

tol = 1e-12;
n = length(T) - length(P);
M = Greville(T,n);
S = LinZero(M,P);
F = DeBoor(P,T,S);

i = find(max(abs(F))> tol);

if isempty(i)
    Z = S;
else
    [P,T] = KnotIns(P,T,S(i));
    Z = Zero(P,T);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function M = Greville(T,n)

c = cumsum(T);
M = (c(n:end-1) - c(1:end-n))/(n-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [P,T] = KnotIns(P,T,S)

n = length(T) - length(P);

for s = S
    T = [T(1)+[-n:-1] T T(end)+[1:n]];
    P = [zeros(1,n) P zeros(1,n)];
    m = length(P);
    k = max(find(T <= s));
   
    j = k-n+2:k;
    
    W = (s - T(j))./(T(j+n-1) - T(j));
    B = W.*P(j) + (1-W).*P(j-1);
   
    P = [P(1:k-n+1) B P(k:m)];
    P = P(n+1:end-n);
    T = [T(n+1:k) s T(k+1:end-n)];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function S = LinZero(M,P)

k = find(P(1:end-1) .* P(2:end) <= 0);

S = M(k) - P(k).*(M(k+1)-M(k))./(P(k+1)-P(k));