Of course, what the Picard iteration is doing in this case is grinding out, one term at a time, the Maclaurin series of the solution. That could be obtained more simply with dsolve(..., type=series).
As for guessgf, how does it come up with an answer? Let's use infolevel[gfun] to get information on what it did.
> infolevel[gfun]:= 2;
guessgf: Trying to find a rational generating function
guessgf: Trying to find an hypergeometric generating function
guessgf: Trying to find an algebraic generating function
guessgf: Trying to find a linear differential equation
listtodiffeq: The ogf seems to satisfy x^2+y(x)-diff(y(x),x)
guessgf: Trying to solve it
[-2 - 2 x - x + 2 exp(x), ogf]
Surprise! It did it by solving the same differential equation we started with.