I know that for the case of linear homogeneous recurrences with constant coefficients, there is a nice method to solve them using generating function. Is there a similar method to solve such recurrences with non-constant coefficients, like the one below?
First note that $p_n$ counts the number of permutations $\sigma$ of $[n]$ such that $\sigma$ has only cycles of length one or two in its cycle decomposition. The recurrence can then be derived by considering whether $1$ is fixed by $\sigma$ or not.
Then note that the recurrence may be written in a more convenient form as $$ p_{n+2}=p_{n+1}+(n+1)p_n\quad (n\ge0)\tag{1} $$ Let $F(x)=\sum_{n\geq0} p_n x^n/n!$. Then $F'(x)=\sum_{n\geq 0} p_{n+1} x^n/n!$ i.e. differentiating corresponds to a left shift of the sequence. Thus translating (1) into EGFs gives us that $$ F''(x)=F'(x)+(xF(x))'\quad F(0)=1, F'(0)=1\tag{2} $$ Thus it suffices to solve the above differential equation. I will leave this task to you.