NB. Coefficients of the Legendre Polynomials (& their derivatives)
NB. Usage: LegP_coeff n --> 2x(n+1) matrix: 
NB. 1st row are coefficients of P_n; the 2nd row is P_n'.
NB. Values are exact rationals: i.e., _315r32
Pshift=: monad : '}: 0,y'
NB. Multiply by x => Shift coefficients to right
acf=: monad : '((2x*y)+1x)%(y+1x)'  NB. (2n+1)/(n+1)
bcf=: monad : 'y %(y+1x)'           NB.      n/(n+1)
LegP_coeff=: monad define
nn=. >: y
if. y=0 do. 
   pc=. 1x                 NB. P_0 = 1
   pcp=. 0                 NB. P_0' = 0
elseif. y=1 do.
   pc=. 0 1x               NB. P_1 = 0 + x
   pcp=. 1x 0              NB. P_1' = 1
elseif. y>1 do.
   pc0=. 1x (0)} pc=. nn $ n=. 0 
   pc1=. 1x (1)} pc 
   while. (n=. >: n) < (nn-1) do.
      pc=. ((acf n)*Pshift pc1) - ((bcf n)*pc0)
      pc0=. pc1
      pc1=. pc
   end. 
end.
pcp=. (}. (i. nn)* pc), 0   NB. coeff of P_n' (dP_n/dx)
(2,nn)$ pc,pcp
)
NB. Evaluate Legendre polynomial: n Pn x --> P_n(x)
Pn =: dyad : '(0{ LegP_coeff x) p. y'
NB. Evaluate derivative of Pn: n dPn x --> d P_n(x)/dx
dPn=: dyad : '(1{ LegP_coeff x) p. y'
NB. For example, try
NB. require 'plot'
NB. load 'local_defs'
NB. x=. _1 + 2*int01 1000
NB. 'pensize 3' plot x;(13 Pn x)

NB. Roots & weights for Gaussian quadrature  
NB. over arbitrary interval [a,b].
NB. Usage: n Gaussian_Quad (a,b) --> roots, wts
Gaussian_Quad=: dyad define
rts=. /:~  >1{ p.  0{ LegP_coeff x
NB. weights: see Abram. & Stegun, p 888.
wts=. 2 % (1- *:rts)* *:(x dPn rts)
NB. transform from [-1,1] to [a,b]:
'a b'=. y
c1=. -: b-a
c2=. -: a+b
rts=. c2 + c1*rts
wts=. c1*wts
|: (2,(x))$ rts,wts
)