function [p, d1, d2] = plate_rr_legendre(order, x)
% PLATE_RR_LEGENDRE  Legendre polynomials and their first two derivatives.
%
%   [p, d1, d2] = plate_rr_legendre(order, x) tabulates P_n(x), P_n'(x) and
%   P_n''(x) for n = 0..order at the points of the row vector x in [-1, 1],
%   as (order + 1) x numel(x) arrays, from the three-term recurrence
%     (n + 1) P_{n+1} = (2 n + 1) x P_n - n P_{n-1}
%   and its derivatives
%     P'_{n+1} = P'_{n-1} + (2 n + 1) P_n,   P''_{n+1} = P''_{n-1} + (2 n + 1) P'_n.

x = x(:)';
m = numel(x);
p = zeros(order + 1, m);
d1 = zeros(order + 1, m);
d2 = zeros(order + 1, m);
p(1, :) = 1;
if order >= 1
    p(2, :) = x;
    d1(2, :) = 1;
end
for n = 1:order - 1
    p(n + 2, :) = ((2 * n + 1) * x .* p(n + 1, :) - n * p(n, :)) / (n + 1);
    d1(n + 2, :) = d1(n, :) + (2 * n + 1) * p(n + 1, :);
    d2(n + 2, :) = d2(n, :) + (2 * n + 1) * d1(n + 1, :);
end
end