<s>[INST] Solving a combintorical problem using generating function 
In how many ways $n$ balls can be divided into groups where each group may contain only one ball or two.  

Each ball can be a singleton, or coupled with one of the $n-1$ options. Hence,
$a_0 = 1$
$a_1 = 1$
$a_n = a_{n-1} + (n-1)a_{n-2}$   
The correspond generating function:  
$$F(x) = 1 + x + \sum\nolimits_{n \ge 2} {\left( {{a_{n - 1}} + (n - 1) \cdot {a_{n - 2}}} \right)}  \cdot {x^n}$$
$$F(x) = 1 + x + \sum\nolimits_{n \ge 2} {{a_{n - 1}}} {x^n} + \sum\nolimits_{n \ge 2} {(n - 1)}  \cdot {a_{n - 2}} \cdot {x^n}$$
$$F(x) = 1 + x + x \cdot \sum\nolimits_{n \ge 1} {{a_n}} {x^n} + {x^2} \cdot \sum\nolimits_{n \ge 0} {(n - 1)}  \cdot {a_n} \cdot {x^n}$$
$$F(x) = 1 + x + x\left( {F(x) - 1} \right) + {x^2}\left( ? \right)$$
How to treat the last expression? $$\sum\nolimits_{n \ge 0} {(n - 1)}  \cdot {a_n} \cdot {x^n}$$
 [/INST] To solve the recurrence that you have given,
$$
\begin{align}
f(x)
&=\sum_{k=0}^\infty a_kx^k\\
&=1+x+\sum_{k=2}^\infty (\color{#C00000}{a_{k-1}}+\color{#00A000}{(k-1)a_{k-2}})x^k\\
&=1+x+\color{#C00000}{x(f(x)-1)}+\color{#00A000}{x^2(xf'(x)+f(x))}\\[10pt]
-1
&=x^3f'(x)+f(x)(x^2+x-1)
\end{align}
$$
This equation can be solved using the integrating factor
$$
\begin{align}
g(x)&=xe^{-\frac1{x^{\vphantom{2}}}+\frac1{2x^2}}\\
g'(x)&=\left(1+\frac1x-\frac1{x^2}\right)e^{-\frac1{x^{\vphantom{2}}}+\frac1{2x^2}}
\end{align}
$$
which yields
$$
-\frac1{x^2}e^{-\frac1{x^{\vphantom{2}}}+\frac1{2x^2}}
=(gf)'(x)
$$
giving
$$
\begin{align}
f(x)
&=\left(\mathrm{erfi}\left(\frac{1-x}{x\sqrt2}\right)+i\right)\frac{\sqrt{2\pi}}{2x}e^{-\frac{(1-x)^2}{2x^2}}\\[6pt]
&=1+x+2x^2+4x^3+10x^4+26x^5+76x^6+232x^7+764x^8+\dots
\end{align}
$$
The series above was verified in Mathematica 8 using
Series[(Erfi[(1-x)/x/Sqrt[2]]+I)Sqrt[Pi/2]/x Exp[-(1-x)^2/2/x^2],{x,0,8},Assumptions->Element[x,Reals]&&x>=0]


To get the exponential generating function,
$$
\begin{align}
g(x)
&=\sum_{k=0}^\infty a_k\frac{x^k}{k!}\\
g'(x)
&=1+\sum_{k=2}^\infty(\color{#C00000}{a_{k-1}}+\color{#00A000}{(k-1)a_{k-2}})\frac{x^{k-1}}{(k-1)!}\\
&=1+\color{#C00000}{g(x)-1}+\color{#00A000}{xg(x)}\\[9pt]
&=(1+x)g(x)
\end{align}
$$
which is easily solved by
$$
\begin{align}
g(x)
&=e^{x+x^2/2}\\
&=1+x+2\frac{x^2}{2!}+4\frac{x^3}{3!}+10\frac{x^4}{4!}+26\frac{x^5}{5!}+76\frac{x^6}{6!}+232\frac{x^7}{7!}+764\frac{x^8}{8!}+\dots
\end{align}
$$
The series above was verified with Mathematica 8 using
CoefficientList[Series[Exp[x+x^2/2],{x, 0, 8}],x]Table[n!,{n, 0, 8}]

</s>