`cheb2leg(::Matrix)` should do tensor product Cheb -> tensor product Leg

Note that chebyshevtransform(::Matrix) is a 2D transform:

julia> using FastTransforms, ClassicalOrthogonalPolynomials

julia> n,m = 3,4; f = (x,y) -> chebyshevt(n, x)chebyshevt(m,y)
#13 (generic function with 1 method)

julia> x = ChebyshevGrid{1}(10);

julia> chebyshevtransform(f.(x', x))
10×10 Matrix{Float64}:
 0.0   1.79736e-18  0.0  1.62102e-17  0.0   9.42055e-18  0.0   1.70874e-18  0.0  -9.27717e-18
 0.0   0.0          0.0  0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0   1.11341e-18  0.0  7.41045e-17  0.0  -4.09803e-19  0.0   1.07597e-17  0.0   5.84104e-19
 0.0   0.0          0.0  0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0   5.90641e-17  0.0  1.0          0.0   1.50729e-16  0.0   2.84217e-16  0.0   8.12948e-17
 0.0   0.0          0.0  0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0  -2.19853e-17  0.0  3.47114e-16  0.0  -1.88411e-17  0.0  -1.22156e-17  0.0  -1.4348e-17
 0.0   0.0          0.0  0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0   6.63185e-18  0.0  2.18604e-16  0.0   1.14231e-17  0.0  -4.61608e-19  0.0   5.85226e-18
 0.0   0.0          0.0  0.0          0.0   0.0          0.0   0.0          0.0   0.0

On the other hand, the following applies 1D transform to each column:

julia> n,m = 3,4; f = (x,y) -> legendrep(n, x)legendrep(m,y)
#15 (generic function with 1 method)

julia> cheb2leg(chebyshevtransform(f.(x', x)))
10×10 Matrix{Float64}:
 0.0  1.83296e-17  0.0   6.46249e-17  0.0   7.18753e-18  0.0   1.30004e-18  0.0  -1.21711e-18
 0.0  0.0          0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0  1.19456e-17  0.0  -8.58383e-17  0.0   9.06258e-19  0.0   2.95814e-17  0.0   4.36288e-19
 0.0  0.0          0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0  0.375        0.0   0.625        0.0   6.39255e-17  0.0   1.35064e-16  0.0   1.36998e-16
 0.0  0.0          0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0  5.4873e-17   0.0   1.12785e-16  0.0  -4.54724e-17  0.0  -1.78334e-17  0.0  -1.84471e-17
 0.0  0.0          0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0  1.45648e-16  0.0   2.70896e-16  0.0   7.99515e-18  0.0   9.33654e-18  0.0   3.97946e-18
 0.0  0.0          0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0

The following gives a quick work-around that I think should be the default behaviour:

julia> function cheb2leg2d(A::AbstractMatrix)
       B = cheb2leg(A); cheb2leg(B')'
       end
cheb2leg2d (generic function with 1 method)

julia> cheb2leg2d(chebyshevtransform(f.(x',x)))
10×10 adjoint(::Matrix{Float64}) with eltype Float64:
 0.0  -2.15114e-17  0.0   9.6853e-17   0.0   1.35234e-17  0.0   4.64116e-18  0.0  -3.28115e-18
 0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0   6.13301e-17  0.0  -1.44889e-16  0.0  -3.0644e-17   0.0   7.0058e-17   0.0   1.17616e-18
 0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0  -2.34741e-17  0.0   1.0          0.0  -5.64224e-17  0.0   1.4927e-16   0.0   3.69325e-16
 0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0  -4.39458e-18  0.0   2.26976e-16  0.0  -6.76889e-17  0.0  -1.92562e-17  0.0  -4.97304e-17
 0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0
 0.0  -1.88093e-17  0.0   4.23769e-16  0.0   4.91028e-18  0.0   1.72571e-17  0.0   1.0728e-17
 0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0          0.0   0.0

Start by locating the cheb2leg(::Matrix) and chebyshevtransform(::Matrix) entry points. Compare the current column-wise behavior with the cheb2leg2d workaround shown in the issue, then add coverage for a tensor-product Chebyshev input. Done means matrix conversion applies in both dimensions and produces the expected tensor-product Legendre coefficients.