I tried to define a simple rule definig how \[Lambda] acts on \[Psi][n]: myrule1 = \[Lambda] \[Psi][n_] -> \[Alpha][n + 1] \[Psi][n + 1]; The result I get is correct provided there's just one \[Lambda] on the RHS of \[Psi][n]. For instance: \[Lambda]^2 \[Psi][n] //. myrule1 isn't computed at all. On the other hand, if I do it step by step: \[Lambda] \[Alpha][1 + n] \[Psi][1 + n] /. myrule1 I get the correct result. I tried to define a new rule: myrule2 = \[Lambda]^m_ \[Psi][n_] -> \[Alpha][n + 1] \[Lambda]^( m - 1) \[Psi][n + 1]; but it doesn't work. Since the recursive method seemed to work, I created a function which multiplies \[Psi][m] by \[Lambda] n times: times\[Lambda][n_] := Nest[Times[\[Lambda], #] /. myrule1 &, \[Psi][m], n] & But this is a very crude way of solving this problem. Do you have any other ideas?