function [x]=levalg(u,b)
[x,~] = levinson(u,b);
end

function [xk,yk] = levinson(u,b)

xk=b(1); yk = u(2);

for k=1:length(u)-2;
    sigma = 1-dot(u(2:k+1),yk);
    
    %computation of x
    t= (b(k+1)-dot(u(k+1:-1:2),xk))/sigma;
    x= [xk-t * yk(k:-1:1); t];
    
    %computation of y, i.e. sol of T(u)y=u
    s= (u(k+2)-dot(u(k+1:-1:2),yk))/sigma;
    y= [yk-s*yk(k:-1:1); s];
    
    xk=x;
    yk=y;
end
    k=length(u)-1;
    sigma = 1-dot(u(2:k+1),yk);
    t= (b(k+1)-dot(u(k+1:-1:2),xk))/sigma;
    xk= [xk-t * yk(k:-1:1); t];
end