$G(n) = 1 + \sum\limits_{ \substack{d|n\\d != 1\\d\neq n}}G(n /d)$, how to solve for $G(n)=n$.

This sequence is A002033 in the OEIS. For whatever reason most of what is below is not in that entry.

Assuming that $G(1) = 1$ the recurrence relation can be rewritten

$$\sum_{d | n} G(d) = \begin{cases} 2 G(n) &\mbox{if } n \ge 2 \\ 1 &\mbox{if } n = 1 \end{cases}.$$

This means the LHS is a Dirichlet convolution. Taking Dirichlet series, if $F(s) = \sum_{n \ge 1} \frac{G(n)}{n^s}$ then we get that

$$F(s) \zeta(s) = 2 F(s) - 1$$

which gives $F(s) (2 - \zeta(s)) = 1$, hence

$$F(s) = \frac{1}{2 - \zeta(s)} = \frac{1}{2} \sum_{k \ge 0} \left( \frac{\zeta(s)}{2} \right)^k.$$

Now, $\zeta(s)^k$ is the Dirichlet series of the function $d_k(n)$, the number of ways to write $n$ as the product of $k$ positive integers. This gives

$$G(n) = \frac{1}{2} \sum_{k \ge 0} \frac{d_k(n)}{2^k}.$$

Somewhat more explicitly, if $n = \prod_i p_i^{e_i}$ is the prime factorization of $n$, then

$$d_k(n) = \prod_i {e_i + k - 1 \choose k-1}$$

and hence

$$\boxed{ G(n) = \frac{1}{2} \sum_{k \ge 0} \frac{1}{2^k} \prod_i {e_i + k - 1 \choose k-1} }.$$

This recovers Sil's claim in the comments that $G$ only depends on the exponents in the prime factorization of $n$. For any particular $n$ the expression $d_k(n)$ is a polynomial in $k$ of degree $\sum_i e_i$ and so a closed form for the above sum is available, basically by leveraging the identity

$$\frac{1}{(1 - x)^{e+1}} = \sum_{k \ge 0} {e + k \choose k} x^k$$

and substituting $x = \frac{1}{2}$ to get

$$2^{e+1} = \sum_{k \ge 0} \frac{1}{2^k} {e + k \choose k}.$$

This result is already enough to compute $G(n)$ in the prime power case and it gives

$$G(p^e) = 2^{e - 1}$$

as also noted by Sil in the comments.

In general the problem reduces to expressing $f(k) = \prod_i {k + e_i -1 \choose e_i}$ as a sum of polynomials of the form ${k + e \choose e}$. There is a straightforward algorithm for doing this for fixed $e_i$ by computing finite differences of the $f(k)$ as follows. Define the backwards difference operator

$$(\Delta f)(k) = f(k) - f(k-1).$$

Observe that $\Delta {k + e \choose e} = {k + e-1 \choose e-1}$. Suppose that $f(k) = \sum_{e=0}^d c_e {k + e \choose e}$ where $d$ is the degree of $f$ as a polynomial in $k$. Then taking backwards differences $r$ times gives

$$(\Delta^r f)(k) = \sum_{e=r}^d c_e {k + e - r \choose e - r}.$$

Substituting $k = 0$ gives

$$(\Delta^r f)(0) = c_r + c_{r+1} + \dots + c_d$$

from which it follows that

$$c_e = (\Delta^{e+1} f)(0) - (\Delta^e f)(0)$$

which gives

$$f(k) = \sum_{e=0}^d \left( (\Delta^{e+1} f)(0) - (\Delta^e f)(0) \right) {k+e \choose e}$$

and hence

$$\sum_{k \ge 0} \frac{f(k)}{2^k} = \sum_{e=0}^d \left( (\Delta^{e+1} f)(0) - (\Delta^e f)(0) \right) 2^{e+1}.$$

---

An alternate approach is the following. We can also write the Dirichlet series as

$$F(s) = \frac{1}{1 - (\zeta(s) - 1)} = \sum_{k \ge 0} (\zeta(s) - 1)^k.$$

$(\zeta(s) - 1)^k$ is the Dirichlet series for a function I'll call $d_k^{+}(n)$, the number of ways to write $n$ as the product of $k$ positive integers each of which is greater than $1$. This gives

$$\boxed{ G(n) = \sum_{k \ge 0} d_k^{+}(n) }.$$

which among other things is 1) a finite sum ($d_k^{+}(n) = 0$ as soon as $k$ is greater than the sum $\sum_i e_i$ of the exponents in the prime factorization of $n$) and 2) clearly an integer. This gives that $G(n)$ is the number of ways to write $n$ as a product of any number of positive integers greater than $1$, or equivalently the number of "Gozinta chains" of $n$ as in Project Euler 548, which is presumably where this question comes from.

We can write $d_k^{+}(n)$ in terms of $d_k(n)$ using inclusion-exclusion: doing so should more or less lead to the first formula.

This expression is also enough to conclude that $G(n)$ only depends on the exponents in the prime factorization of $n$, and actually it suggests the following alternate generating function expression for $G$. Write

$$G \left( \prod_{i=1}^m p_i^{e_i} \right) = g(e_1, e_2, \dots e_m).$$

Then $g$ has multivariate generating function

$$\sum_{e_i \ge 0} g(e_1, e_2, \dots e_m) x_1^{e_1} x_2^{e_2} \dots x_m^{e_m} = \frac{1}{2 - (\sum_{e_i \ge 0} x_1^{e_1} x_2^{e_2} \dots x_m^{e_m})} = \frac{1}{2 - \prod_{i=1}^m \frac{1}{1 - x_i}}$$

so this multivariate generating function is even rational: it can be rewritten

$$\displaystyle \frac{\prod_{i=1}^m (1 - x_i)}{2 \prod_{i=1}^m (1 - x_i) - 1}.$$

Written this way the calculations for small values of $m$ are more transparent. For $m = 1$ we get

$$\sum_{e \ge 0} g(e) x^e = \frac{1 - x}{1 - 2x}$$

recovering $g(e) = 2^{e-1}$, and for $m = 2$ we get

$$\sum_{e_1, e_2 \ge 0} g(e_1, e_2) x_1^{e_1} x_2^{e_2} = \frac{1}{2 - \frac{1}{(1 - x_1)(1 - x_2)}} = \frac{(1 - x_1)(1 - x_2)}{1 - 2x_1 - 2x_2 + 2x_1 x_2}$$

from which both the computations of $g(1, e)$ and $g(e, e)$ recorded in Sil's community wiki answer can be recovered. To get a sense of what can be done here, let's for a moment consider the generating function above as a generating function in $x_1$ only, with coefficients which are generating functions in $x_2$. This gives

$$\frac{1 - x_1}{2 - 2x_1 - \frac{1}{1 - x_2}} = \frac{1 - x_1}{\frac{1 - 2x_2}{1 - x_2} - 2x_1}.$$

Factoring out $\frac{1 - 2x_2}{1 - x_2}$ from the denominator gives

$$\frac{1 - x_2}{1 - 2x_2} \frac{1 - x_1}{1 - \frac{2 - 2x_2}{1 - 2x_2} x_1} = (1 - x_1) \sum_{k \ge 0} 2^k \frac{(1 - x_2)^{k+1}}{(1 - 2x_2)^{k+1}} x_1^k.$$

Extracting the coefficient of $k$ then gives, for $k \ge 1$,

$$\sum_{e \ge 0} g(k, e) x^e = 2^k \frac{(1 - x)^{k+1}}{(1 - 2x)^{k+1}} - 2^{k-1} \frac{(1 - x)^k}{(1 - 2x)^k}$$

or, simplifying slightly,

$$\boxed{ \sum_{e \ge 0} g(k, e) x^e = 2^{k-1} \frac{(1 - x)^k}{(1 - 2x)^{k+1}} }.$$

For $k = 1$ this recovers $g(1, e) = (e + 2) 2^{e-1}$ as in Sil's answer. For general $k$ this gives

$$\boxed{ g(k, e) = 2^{k-1} \sum_{i=0}^k (-1)^i 2^{e-i} {k \choose i} {k + e - i \choose k} }.$$