function P = legendre_p(n)
  % Return the `n`th Legendre polynomial.
  %
  % INPUT:
  %
  %   * ``n`` - The index of the polynomial that we want.
  %
  % OUTPUT:
  %
  %   * ``P`` - A polynomial function of one argument.
  %
  if (n < 0)
    % Can't do anything here. Return nothing.
    P = NA;
  elseif (n == 0)
    % One of our base cases.
    P = @(x) 1;
  elseif (n == 1)
    % The second base case.
    P = @(x) x;
  else
    % Not one of the base cases, so use the recursive formula.
    prev = legendre_p(n-1);
    prev_prev = legendre_p(n-2);
    P = @(x) (1/n).*( (2*n - 1).*x.*prev(x) - (n-1).*prev_prev(x) );
  end
end
