function [Q,N] = nubs(CP, knots, p)

n = length(knots) - 5;
u = linspace(knots(4), knots(n+2), p);
d = size(CP,2);

Q1 = zeros(n-1,d,p);

N = blends(knots,u);  %to be defined

for j = 1:p
    
    for i = 2:n-1
        
        %Q1(i,:,j) = zeros(1,d);
    
        for k = 0:3
        
            for l = 1:d
                Q1(i,l,j) = Q1(i,l,j) + N(i+k-1,4,j)*CP(i+k-1,l);
            end
        end
    end
end

Q = Q1(2:n-1,:,:);
%plotting commands