Single-sided Z transform with difference equations and the system function I am working on this problem:
Given an impulse response,

find the system function,
find the difference equation representation,
find pole-zero plot,
find output $y[n]$ if the input is $x[n] = 0.25^n * u[n]$

Here is what I have so far (MATLAB code)

Find H(z)
$$H(z) = \frac{5}{1 - 0.25z^{-1}}\quad,\qquad|z| &gt; 0.25$$

For part 2 I used the fact that
$$H(z)=\frac{Y(z)}{X(z)}$$
I cross multiplied, then took the inverse $\mathcal Z$ transform and got:
$$y[n] - 0.25y[n] = 5x[n]$$

For part 3 I used the zplane function:
 
 b = [5 0];
 a = [1 -0.25];
 figure; zplane(b,a);

Part 4 is where I am getting confused. I took the one-sided $\mathcal Z$-transform of the difference equation to get
$$Y^+(z)-0.25[y[-1]+z^{-1}Y^+(z)] = 5X^+(z)$$
Assuming the system is casual (due to unit-step function in impulse response and input $x[n]$, also no initial conditions given in the problem),
$$y[-1] = 0$$
So, after rearranging, $$Y^+(z)=\frac{5}{1-0.5z^{-1}+0.0625z^{-2}}\quad, \qquad |z| &gt; 0.25$$
Using the residuez function for partial fraction decomp
 
 b1 = [5 0 0];
 a1 = [1 -0.5 0.0625];
 [R p C] = residuez(b1,a1)
 
 R =

      0
      5


 p =

     0.2500
     0.2500


 C =

      0



Then taking the inverse $\mathcal Z$-transform of the terms generated by residuez,
$$y[n] = 5(1/4)^nu[n]$$
However, when I check, this is not the case. For the check, I first create a MATLAB generated $y[n]$ sequence using filter
 
%MATLAB check for part 4
n = 0:50;  % compare first 50 samples
x = (1/4).^n;  
y = filter(b,a,x);
Then, create my check sequence calculated from the inverse $\mathcal Z$-transform and the residuez function

ycheck =  5*(0.25).^n;
error = max(abs(y-ycheck))
The sequences are not the same, and error nowhere close to nominal. Can anyone tell me where I went wrong please?
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.