If you have a difference equation like this:

    eqn = y[n + 2] - 5 y[n + 1] + 6 y[n] == 
      3 (0.5)^(1 + n) UnitStep[1 + n] + 5 (0.5)^n UnitStep[n]

6 y[n] - 5 y[1 + n] + y[2 + n] == 
 5 0.5^n UnitStep[n] + 3 0.5^(1 + n) UnitStep[1 + n]

With the boundary conditions:

y[-1] = 11/6;

y[-2] = 37/36;

You will find that it is easily and correctly solvable using RSolve command. However, if you try to solve it using z-transforms, the process is less clear and my methodology does not yield the correct answer. First take the z-transform of eqn:

ztde = ZTransform[eqn, n, z]

-z^2 y[0] - z y[1] + 6 ZTransform[y[n], n, z] + z^2 ZTransform[y[n], n, z] - 
  5 (-z y[0] + z ZTransform[y[n], n, z]) == (13. z)/(-1. + 2. z)

As can be seen the expression ztde of the z-transformed version of eqn is in terms of boundary conditions y[1] and y[0]. These are not the boundary conditions supplied. Examining this problem I thought that the way to get these boundary conditions from the y[-1] and y[-2] supplied would be to use the zero input response part (some older differential equation books would call this the 'homogeneous part') of eqn, use n=-1 and n=-2, and solve for y[0] and y[1]. That would be

boundaryvaleqn = eqn[[1]] == 0

6 y[n] - 5 y[1 + n] + y[2 + n] == 0

boundaryval1 = boundaryvaleqn /. n -> -1

11 - 5 y[0] + y[1] == 0

boundaryval2 = boundaryvaleqn /. n -> -2

-3 + y[0] == 0

newboundval = Solve[{boundaryval1, boundaryval2}, {y[1], y[0]}]

{{y[1] -> 4, y[0] -> 3}}

Applying new boundary values to ztde

zeqn = ztde /. newboundval

{-4 z - 3 z^2 + 6 ZTransform[y[n], n, z] + z^2 ZTransform[y[n], n, z] - 
   5 (-3 z + z ZTransform[y[n], n, z]) == (13. z)/(-1. + 2. z)}

Solving for the z-transform and inverting

zsolve = Solve[zeqn, ZTransform[y[n], n, z]]

{{ZTransform[y[n], n, z] -> (-11. z + 3. z^2 + (13. z)/(-1. + 2. z))/(
   6. - 5. z + z^2)}}

solution = y[n] /. InverseZTransform[zsolve, z, n][[1]]

0.0666667 2.^(-1. n) (26. + 5. 2.^(1. + 2. n) + 2.^n 3.^(2. + n))

This is an answer but it is not the correct answer. My guess is that something is wrong in my assumptions in recalculating the boundary conditions or is there some way so when the z-transform is taken there is a way to specify the correct boundary conditions. Any thoughts or ideas in approaching this?