I think you didn't read all the way to the end of the doc for ORPOL. It includes an example of how to generate Legendre polynomials by using the usual three-term recurrence. The doc example computes the Legendre polynomials up through degree 6. The only thing I added was writing to a SAS data set and plotting the polynomials on [-1, 1]:
proc iml;
maxDegree = 6;
/* evaluate polynomials at these points */
x = T( do(-1,1,0.05) );
/* define the standard Legendre Polynomials
Using the 3-term recurrence with
A[j]=0, B[j]=(2j-1)/j, and C[j]=(j-1)/j
and the standardization P_j(1)=1
which implies P_0(x)=1, P_1(x)=x. */
legendre = j(nrow(x), maxDegree+1);
legendre[,1] = 1; /* P_0 */
legendre[,2] = x; /* P_1 */
do j = 2 to maxDegree;
legendre[,j+1] = (2*j-1)/j # x # legendre[,j] -
(j-1)/j # legendre[,j-1];
end;
*print legendre;
L = x || Legendre;
create Legendre from L[c=('x' || ('L0':'L6'))];
append from L;
close;
QUIT;
proc sgplot data=Legendre;
series x=x y=L0;
series x=x y=L1;
series x=x y=L2;
series x=x y=L3;
series x=x y=L4;
series x=x y=L5;
series x=x y=L6;
run;