function [t,p] = lagrange(x,y,a,b);
%
% Poynomial interpolation in Lagrange Form. 
%
% Input:   x : vector with the nodes
%          y : Vector with the corresponding y-values. 
%      [a,b] : Interpolation interval 
% Output:  t : vector of equally spaced points in [a,b]
%          p : the value of the interpolation 
%              polynomial in [a,b].

n = length(x);
t = linspace(a,b,1001)';

p = zeros(size(t));
L = cardinal(x,t);
for k=1:n
    p = p + y(k)*L(:,k);
end


function L = cardinal(x,t)
%
% The cardinal functions lk of x, evaluated in t. 
% t is a column vector. 
% Output: L = [l1(t),l2(t),....,ln(t)]
% 
n = length(x);
L = ones(length(t),n);
for k=1:n
    for i=1:n
        if i ~= k
            L(:,k) = L(:,k).*(t-x(i))/(x(k)-x(i));
        end
    end
end
