Help with Lagrange Interpolation: How can I graph the polynomial of my data sett?

I am attempting to use Lagrange interpolation to find the polynomial for my data set. How can I alter my code to do this? Below is the code so far, but it only plots the points:

x = [1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010];
y = [76.21 92.23 106.0 123.2 132.2 151.3 179.3 203.3 226.5 248.8 281.4 308.7];

disp([x; y])
u = 1880:2020:12.25;
v = polyinterp(x,y,u);
plot(x,y,'o',u,v,'r-');

function v = polyinterp(x,y,u)
    n = length(x);
    v = zeros(size(u));
    for k = 1:n
    w = ones(size(u));
    for j = [1:k-1 k+1:n]
    w = (u-x(j))./(x(k)-x(j)).*w;
    end
    v = v + w*y(k);
    end
end