3+ dim calculus operator do not work and Syntax proposal

Currently,

julia> Derivative(Chebyshev()^2, [1,0])
DerivativeWrapper : Chebyshev() ⊗ Chebyshev() → Ultraspherical(1) ⊗ Chebyshev()
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  2.0  0.0  0.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  2.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  3.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
  ⋮    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱   ⋱

julia> Integral(Chebyshev()^2, [1,0])
ERROR: Implement Integral(Chebyshev() ⊗ Chebyshev(),[1, 0])
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] DefaultIntegral(sp::TensorSpace{Tuple{Chebyshev{ChebyshevInterval{Float64}, Float64}, Chebyshev{ChebyshevInterval{Float64}, Float64}}, DomainSets.FixedIntervalProduct{2, Float64, ChebyshevInterval{Float64}}, Float64}, k::Vector{Int64})
   @ ApproxFunBase ~/.julia/packages/ApproxFunBase/VD9Rj/src/Operators/banded/CalculusOperator.jl:41
 [3] Integral(::TensorSpace{Tuple{Chebyshev{ChebyshevInterval{Float64}, Float64}, Chebyshev{ChebyshevInterval{Float64}, Float64}}, DomainSets.FixedIntervalProduct{2, Float64, ChebyshevInterval{Float64}}, Float64}, ::Vararg{Any})
   @ ApproxFunBase ~/.julia/packages/ApproxFunBase/VD9Rj/src/Operators/banded/CalculusOperator.jl:53
 [4] top-level scope
   @ REPL[94]:1

Integral(Chebyshev()) ⊗ I
KroneckerOperator : Chebyshev() ⊗ ApproxFunBase.UnsetSpace() → Chebyshev() ⊗ ApproxFunBase.UnsetSpace()

It would be nice to have a Syntax for Integral(Chebyshev()^2, [1,0]) for integration over the first variable as well as for DefiniteIntegral.

Furthermore, currently the following throws an error:

julia> Op = Integral(Chebyshev()) ⊗ I ⊗ I 
KroneckerOperator : Chebyshev() ⊗ ApproxFunBase.UnsetSpace() ⊗ ApproxFunBase.UnsetSpace() → Chebyshev() ⊗ ApproxFunBase.UnsetSpace() ⊗ ApproxFunBase.UnsetSpace()

julia> f3 = Fun(Chebyshev()^3, rand(20));

julia> Op * f3
ERROR: Implement Conversion from Chebyshev{ChebyshevInterval{Float64}, Float64} to TensorSpace{Tuple{Chebyshev{ChebyshevInterval{Float64}, Float64}, ApproxFunBase.UnsetSpace}, DomainSets.VcatDomain{2, Float64, (1, 1), Tuple{ChebyshevInterval{Float64}, DomainSets.FullSpace{Float64}}}, Union{}}
Stacktrace:

as well as the same for the Derivative.

Maybe I am misusing the package for something it is not intended to be, but I want to use it for implementing closed form integration, marginalization, and derivatives of SOS polynomials and therefore would need this operations for higher dimensions.

@jishnub Do you have an advice for me where I have to start for implementing this?

Start with the Integral and Derivative entry points and the DefaultIntegral implementation in Operators/banded/CalculusOperator.jl, which appears in the stack trace. Reproduce the tensor-space and three-dimensional examples, then trace the conversion and Kronecker-operator application failures. Done should include higher-dimensional integration and differentiation, including the requested variable-selection syntax and matching behavior for DefiniteIntegral.