VariationOfParameters[homogeneouseqn_, forcingterm_, y_, n_] := Block[{sol, r, y1, y2, c, u1, u2}, sol = RSolve[homogeneouseqn, y, n]; y1[r_] := (y[r] /. sol[[1]] /. {C[1] -> 1, C[2] -> 0}); y2[r_] := (y[r] /. sol[[1]] /. {C[1] -> 0, C[2] -> 1}); c[r_] := Casoratian[{y1[r], y2[r]}, r]; u1 = -Sum[y2[n + 1] forcingterm/c[n + 1], {n, 0, n - 1}]; u2 = Sum[y1[n + 1] forcingterm/c[n + 1], {n, 0, n - 1}]; Simplify[y1[n] u1 + y2[n] u2, Element[n, Integers]] ]
VariationOfParameters[y[n + 2] - 4 y[n] == 0, n^2, y, n]