How to introduce an interpolation function into a neville's algorithm to solve a polynomial interpolation

The question is:

interpolate f: [ 1, 1], f (x) = 1/(25x^2 + 1)

with a Lagrange polynomial at the equidistant points x_k = 1 + 2k/n, k = 0,...,n and plot the graph of pn for n = 10, 20, 40.

This is what I have so far

% plot out the n=10, 20, 40 polynomials which interpolate the 
%    function f(x) = 1/(25x^2+1) on the interval [-1,1] using 
%    equally spaced points
% function [yy]= nev[1/(25*xx.^2+1)]; 

xx = -1:0.1:1;  % plot points in each case

% first plot the original function
yy = 1 ./ (25 * xx.^2 + 1); 
plot(xx,yy);
hold on;               % this makes the plots that follow pile up

for n = [10 20 40]
   x = -1:2/n:1;         
   Q = 1 ./ (25 * x.^2 + 1);   % get the interpolation points

     for i = 1: 201
        yy(i) = nev(xx(i), n, x, Q);  %interpolate
     end

     plot(xx,yy);
  end

hold off

But the function is not interpolating. What should the code be to get the desired graph?