Already done, including differentiation:
julia> P = Jacobi(0.1,0.2)
Jacobi{Float64}
julia> @time P[range(-1,1;length=1000),1:1000] # 1000x1000 Vandermonde matrix
0.009533 seconds (4.03 k allocations: 7.723 MiB)
1000×1000 Array{Float64,2}:
1.0 -1.2 1.32 -1.408 1.4784 -1.53754 1.58879 -1.63418 1.67504 -1.71226 … -4.33032 4.33119 -4.33206 4.33293 -4.3338 4.33467 -4.33554
1.0 -1.1977 1.31274 -1.39289 1.45237 -1.4974 1.53128 -1.55598 1.57277 -1.58255 -0.100404 0.111146 -0.121432 0.131222 -0.140477 0.149161 -0.157239
1.0 -1.1954 1.30549 -1.37784 1.42656 -1.4578 1.47486 -1.47976 1.47385 -1.45815 -0.156531 0.157461 -0.157131 0.155542 -0.152711 0.14866 -0.143423
1.0 -1.19309 1.29826 -1.36287 1.40097 -1.41873 1.41952 -1.40548 1.37819 -1.3389 -0.0688575 0.0554815 -0.0414527 0.0269404 -0.0121192 -0.00283266 0.0177356
1.0 -1.19079 1.29104 -1.34798 1.37561 -1.38018 1.36523 -1.33312 1.28573 -1.22467 -0.074041 0.0859243 -0.0964202 0.105362 -0.112608 0.118045 -0.121586
1.0 -1.18849 1.28384 -1.33315 1.35047 -1.34216 1.31198 -1.26265 1.1964 -1.11533 … -0.0216154 0.00553179 0.0106463 -0.0265952 0.0419959 -0.0565408 0.0699398
⋮ ⋮ ⋱ ⋮
1.0 1.09079 1.12737 1.13802 1.13063 1.10867 1.07422 1.02876 0.97354 0.909709 … 0.067262 0.074961 0.0814524 0.0866333 0.0904222 0.0927598 0.0936102
1.0 1.09309 1.13425 1.15179 1.15349 1.14272 1.12136 1.09064 1.05154 1.00487 0.0370578 0.0263838 0.0154038 0.00424984 -0.00694394 -0.0180431 -0.0289146
1.0 1.0954 1.14115 1.16562 1.17656 1.17726 1.16946 1.15424 1.13237 1.10444 0.115428 0.11449 0.112636 0.109883 0.106253 0.101777 0.0964897
1.0 1.0977 1.14807 1.17952 1.19984 1.21228 1.21854 1.21959 1.21611 1.20853 0.0894469 0.0962089 0.102579 0.108532 0.114045 0.119096 0.123666
1.0 1.1 1.155 1.1935 1.22334 1.2478 1.2686 1.28672 1.30281 1.31728 2.09594 2.09615 2.09636 2.09657 2.09678 2.09699 2.0972
julia> D = Derivative(axes(P,1));
julia> (D*P)[range(-1,1; length=1000),1:1000] # 1000x1000 derivative collocation matrix
(1000×999 view(::Jacobi{Float64}, -1.0:0.002002002002002002:1.0, 1:999) with eltype Float64) * (999×1000 view(::BandedMatrices.BandedMatrix{Float64,Adjoint{Float64,InfiniteArrays.InfStepRange{Float64,Float64}},InfiniteArrays.OneToInf{Int64}}, 1:999, 1:1000) with eltype Float64):
0.0 1.15 -3.63 7.568 -13.0592 20.1802 -28.9954 39.5608 -51.9261 66.136 -82.2312 … 1.78935e6 -1.79331e6 1.79727e6 -1.80124e6 1.80521e6
0.0 1.15 -3.6229 7.53154 -12.9472 19.9131 -28.451 38.5645 -50.2406 63.4521 -78.1577 -2492.57 2366.12 -2230.06 2084.93 -1931.3
0.0 1.15 -3.6158 7.49514 -12.8356 19.6481 -27.913 37.5841 -48.5907 60.8401 -74.2189 123.465 -279.471 433.396 -584.002 730.08
0.0 1.15 -3.60869 7.45882 -12.7246 19.3852 -27.3813 36.6196 -46.9759 58.2986 -70.4114 1187.1 -1221.03 1240.32 -1244.72 1234.16
0.0 1.15 -3.60159 7.42257 -12.6141 19.1243 -26.8558 35.6708 -45.3957 55.8262 -66.7321 -605.409 504.391 -395.189 279.546 -159.308
0.0 1.15 -3.59449 7.3864 -12.5041 18.8655 -26.3365 34.7374 -43.8495 53.4217 -63.178 … 802.793 -783.992 749.469 -699.899 636.262
⋮ ⋮ ⋮ ⋱ ⋮
0.0 1.15 3.43659 6.85739 11.3677 16.8875 23.3057 30.4844 38.2618 46.4568 54.8727 … -363.761 -279.673 -191.019 -99.2159 -5.73067
0.0 1.15 3.44369 6.89249 11.4722 17.1292 23.7846 31.3366 39.663 48.6237 58.0633 920.209 929.935 928.495 915.895 892.276
0.0 1.15 3.4508 6.92767 11.5771 17.3729 24.2694 32.2034 41.0958 50.8533 61.369 288.345 399.365 507.296 611.27 710.452
0.0 1.15 3.4579 6.96293 11.6826 17.6186 24.7601 33.0849 42.5609 53.1468 64.7926 -1531.18 -1425.78 -1314.56 -1197.96 -1076.45
0.0 1.15 3.465 6.99825 11.7885 17.8663 25.2567 33.9812 44.0586 55.5055 68.3371 9.44618e5 9.46611e5 948607.0 9.50605e5 952605.0
(Don’t ask why the collocation matrix is left lazy… it can be materialised but takes 4s right now, could be sped up to roughly 0.01s pretty easily though.)
Thank You again Sheehan