For simple case where the 2 sequences are assumed to start at n=0, then this can be done in the time domain using the ListConvolve in Mathematica and using conv in Matlab.
The harder case is when the origin is not located at the first element in the sequence. I.e. one or both of the sequences is not causal.
Mathematica
case 1
Convolve 2 sequences where the origin is assumed to be at the first element in each sequence.
x1={1,2,3,4,5}; x2={4,5,6}; ListConvolve[x2,x1,{1,-1},0]
{4,13,28,43,58,49,30}
case 2
Convolve 2 sequences where the origin is located at different location from the first element of the sequence, use DiracDelta function, and the DTFT approach.
Example convolve \(x=\{1,2,0,2,1,4,01,2,2\}\) with \(h=\{1/4,1/2,1/4\}\) where the origin of \(h\) is under the second element 1/2.
(*Write down the h and take its DTFT *) Clear[n,w] h = 1/4DiscreteDelta[n+1] + 1/2DiscreteDelta[n]+1/4 DiscreteDelta[n-1]; hTransformed = FourierSequenceTransform[h,n,w]
\[ \frac {1}{4} e^{-i w} \left (1+e^{i w}\right )^2 \]
Write down the x sequence and take its DTFT
x = 1 DiscreteDelta[n]+2 DiscreteDelta[n-1]- 2 DiscreteDelta[n-2]+DiscreteDelta[n-3]+ 4 DiscreteDelta[n-4]-DiscreteDelta[n-5]+ DiscreteDelta[n-6]+2 DiscreteDelta[n-7]; xTransformed = FourierSequenceTransform[x,n,w]
\[ e^{-7 i w} \left (e^{i w}-e^{2 i w}+4 e^{3 i w}+e^{4 i w}-2 e^{5 i w}+2 e^{6 i w}+e^{7 i w}+2\right ) \]
Now multiply the DTFT’s and take the inverse
z = InverseFourierSequenceTransform[ xTransformed hTransformed,w,n]
Now convolve z with h again, where z is the convolution of x and h found above. This can be done as follows in one command
N@InverseFourierSequenceTransform[ xTransformed hTransformed hTransformed,w,n]