function [Nshape,sp]=Bsp(U,p,x);
%
m = length(U)-1;
n = m-p-1;
mm = m+1;
pm = p+1;
%
%Computing the span
if x>=U(n+1)-(U(n+1)*0.0001);
    sp = n+1;
else
    b=x>=U;
    b(b==0)=[];
    sp=(length(b));
end
%
%Calculating the basis functions
%
Nshape = zeros(1,pm);
left = zeros(1,p);
right = zeros(1,p);
Nshape(1,1) = 1;
%
for j = 1:p
    left(j) = (x-U(sp+1-j));
    right(j) = (U(sp+j)-x);
    saved =0;
    r=0;
    while r<j;
        temp = Nshape(1,r+1)/(right(1,r+1)+left(1,j-r));
        Nshape(1,r+1) = saved+right(1,r+1)*temp;
        saved = left(1,j-r)*temp;
        r=r+1;
    end
    Nshape(1,j+1) = saved;
end
