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]