function lkttau = lk( xx, k, ttau)
%
% Lagrange basis function lk(x)
% on the nodes xx
% evaluated at x=ttau (a scalar or a vector of points)
% 
n=size(xx,2)-1;
%
N=size(ttau,2)-1;
lkttau=zeros(1,N+1);
for l=1:N+1
    tau=ttau(1,l);
vv=[xx(1,1:k-1) xx(1,k+1:n+1)];
vv1=tau*ones(1,n)-vv;
vv2=xx(1,k)*ones(1,n)-vv;
lktau=prod(vv1)/prod(vv2);
lkttau(1,l)=lktau;
end
%
end