%function out=LaGrangeDer(xdata,ydata,x)
%
% Input:
%              xdata: points x0...xn
%              ydata: function values f(x0)...f(xn)
%              these should be row vectors
%              x: point at which to interpolate the derivative of the polynomial
% Output:
%              out: This is P'(x), where P is the unique nth-order
%              polynomial with P(xi)=f(xi)

function out=LaGrangeDer(xdata,ydata,x)
n=size(xdata,2)-1;
out=0;

for i=0:n
    out=out+DLP(i,xdata,x)*ydata(i+1);
    %this implements lagrange interpolation: P'(x)=sum_i(L'_i(x)*f(x_i))
end

end

%function out=DLP(i,xdata,x)
%
% Input:
%              i: index number of the LaGrange polynomial to compute
%              xdata: coordinate of the points x0...xn
%              x: point at which to compute L'_i
% Output:
%              out: value of L_i(x)

function out=DLP(i,xdata,x)
n=size(xdata,2)-1;
out=0;
for j=0:n
%   we use generalized product rule, here taking
%   derivative of the jth term of L_i
    if not(eq(i,j))
       DLPterm=1/(xdata(i+1)-xdata(j+1));
       %derivative of (x-x_j)/(x_i-x_j) is 1/(x_i-x_j)
       for k=0:n
          if and(not(eq(i,k)),not(eq(j,k)))
              DLPterm=DLPterm*(x-xdata(k+1))/(xdata(i+1)-xdata(k+1));
              %all the terms except the jth term still appear in the
              %product
          end
       end
       out=out+DLPterm;
    end
end
end