Single-sided Z transform with difference equations and the system function
You misinterpreted the result of residuez. Note that there is also a fourth output argument, namely the multiplicity of the poles (which is $2$ in your example). In your case, the residuez function doesn't help much, because you know the poles already, and the result is already in a partial fractions form:
$$Y(z)=X(z)H(z)=\frac{5}{\left(1-\frac14 z^{-1}\right)^2}\tag{1}$$
You can either look up the inverse $\mathcal{Z}$-transform of $(1)$ in a table, or, maybe even more straightforward, compute the result by convolution in the time domain:
$$y[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]\tag{2}$$
Solving $(2)$ for the given sequences is very straightforward.