I have recently written a code where I use Clenshaw's summation formula with Chebyshev polynomials
$S(x)=\sum_{k=0}^nc_kT_k(x)=b_0+xb_1$
$T_{k+1}(x)=2xT_k(x)-T_{k-1}(x)$
$T_0(x)=1~~~ T_1(x)=x$
$b_{n+2}=b_{n+1}=0$
$b_j=c_j+2xb_{j+1}-b_{j+2};~~~ j = n,n-1,...,1$
I would like to use Clenshaw's algorithm for a two-variable sum S(x,y)
$S(x,y)=\sum_{k=0}^n\sum_{l=0}^mc_{kl}T_k(x)T_l(y)$
I have tried with brute force to just expand the current method as much as possible, but since my understanding of the issue is at minimum, I have not been able to solve it. Any help would be amazing!