% function to interpolate and downgrid the derivatives matrices
% and the identy

function P = Poly_In(xnew,x)

%% Find the weights
np=length(x);
a=min(x); b=max(x);
w=ones(np,1);
C=4/(b-a);
for j=1:np
    for k=1:j-1
        w(j)=w(j)*(x(j)-x(k))*C;
    end
    for k=j+1:np
        w(j)=w(j)*(x(j)-x(k))*C;
    end
end
w=1./w;
%%





P = bsxfun(@minus,xnew,x');          %All xnew(j) - x(k)
P = bsxfun(@rdivide,w',P);           %All w(k)/(xnew(j) - x(k))
P = bsxfun(@rdivide, P, sum(P,2));   % Normalization
P(isnan(P)) = 1;                     % Remove NaNs
end


%% Example
%  [x,~,w] = chebpts(n+2);  % Second order derivative
%  xnew = chebpts(n);       % Smaller grid
%  P = Poly_In(xnew,x);     % Matrix to cast y to ynew by ynew = P*y