f(x) = F[i](x) for i-1 <= x <= i, i = 1 .. infinity, is a solution on [1, infinity). 
As I said, it's trickier to get a solution on the whole real line.  But it will work if
you take F[1](x) to be a smooth function that is 0 in, say, [0, 1/3] union [2/3, 1].
Again, for n > 1, F[n](x) is the solution of y'(x) = y(x) - F[n-1](x-1) with y(n-1) = F[n-1](n-1),
while for n < 1, F[n](x) = F[n+1](x+1)-F[n+1]'(x+1), and you take f(x) = F[i](x) for i-1 <= x <= i.

The real challenge is to get an analytic solution.

f(x) = F[i](x) for i-1 <= x <= i, i = 1 .. infinity, is a solution on [1, infinity). 
As I said, it's trickier to get a solution on the whole real line.  But it will work if
you take F[1](x) to be a smooth function that is 0 in, say, [0, 1/3] union [2/3, 1].
Again, for n > 1, F[n](x) is the solution of y'(x) = y(x) - F[n-1](x-1) with y(n-1) = F[n-1](n-1),
while for n < 1, F[n](x) = F[n+1](x+1)-F[n+1]'(x+1), and you take f(x) = F[i](x) for i-1 <= x <= i.

The real challenge is to get an analytic solution.

No, simplify does not care about "small" or "large".  The imaginary part in this case is exactly 0, because 1/6*120^(1/3)*15^(1/6) = 1/180*120^(2/3)*15^(5/6) = 1/3*15^(1/2).  The nonzero result you get by evaluating it in floating point is entirely due to roundoff error.

If the limit exists, that should give the correct result.  Of course there's also the case x=0 to consider. 
But in general it won't always work: there are examples where the limit along each straight line through 0 exists (and is always the same), but the limit as x and y both go to 0 does not.  Try e.g.

abs(x^2 - y)^(x^2 + y^2)

If the limit exists, that should give the correct result.  Of course there's also the case x=0 to consider. 
But in general it won't always work: there are examples where the limit along each straight line through 0 exists (and is always the same), but the limit as x and y both go to 0 does not.  Try e.g.

abs(x^2 - y)^(x^2 + y^2)

In fact you can let f(x) be any smooth function on an interval of length 1, say [0,1], with f'(1) = f(1) - f(0), and get a solution on [0,infinity), obtaining it iteratively on each interval [n,n+1]. 
For example, with f(x) = x^3 - 2*x^2 on [0,1], I get 3*x-2*x^2+x^3-3*exp(x-1) on [1,2], x^3-2*x^2+6*x+3*x*exp(x-2)-12*exp(x-2)-3*exp(x-1) on [2,3],

> F[1]:= x -> x^3 - 2*x^2;
  for i from 2 to 5 do
      F[i]:= unapply(simplify(rhs(dsolve({D(f)(x) = f(x) - F[i-1](x-1), f(i-1) = F[i-1](i-1)}))), x);
   end do;

It's trickier to get a solution on the whole real line (the backward iteration loses differentiability at the endpoints of the intervals).

In fact you can let f(x) be any smooth function on an interval of length 1, say [0,1], with f'(1) = f(1) - f(0), and get a solution on [0,infinity), obtaining it iteratively on each interval [n,n+1]. 
For example, with f(x) = x^3 - 2*x^2 on [0,1], I get 3*x-2*x^2+x^3-3*exp(x-1) on [1,2], x^3-2*x^2+6*x+3*x*exp(x-2)-12*exp(x-2)-3*exp(x-1) on [2,3],

> F[1]:= x -> x^3 - 2*x^2;
  for i from 2 to 5 do
      F[i]:= unapply(simplify(rhs(dsolve({D(f)(x) = f(x) - F[i-1](x-1), f(i-1) = F[i-1](i-1)}))), x);
   end do;

It's trickier to get a solution on the whole real line (the backward iteration loses differentiability at the endpoints of the intervals).

plots:-implicitplot(eq2, x=-1..1, y=-1..1,gridrefine=3);

plots:-implicitplot(eq2, x=-1..1, y=-1..1,gridrefine=3);

The gfun package was originally developed by Bruno Salvy and Paul Zimmerman.   I don't know who's in charge of it now.  I think you should report it, to make sure that it does get entered into the bug database and (hopefully) eventually fixed.

1) It's a bug

2)

> with(gfun):
    P:= -(5+8*x-x^2)/(-1+3*x^2);
    de:= holexprtodiffeq(P, y(x));

de := (-1+3*x^2)*y(x)+5+8*x-x^2

> diffeqtorec(de, y(x), u(t));

{3*u(t)-u(t+2), u(0) = 5, u(1) = 8, u(2) = 14}

Let c = 0.0001.  For "small" n, S[n] is a good approximation to Int(sin(c*x^2),x = 0 .. n) = 1/2*2^(1/2)*(Pi/c)^(1/2)*FresnelS(n*2^(1/2)*(1/Pi*c)^(1/2))

That explains the rapid rise and the look of the graph for the first few thousand points: note that as n -> infinity this is
1/4*2^(1/2)*Pi^(1/2)/c^(1/2)-1/2*1/c*cos(n^2*c)/n+O(1/(n^3)).

Now exp(I*c*(x+31416)^2) = exp(I*c*31416^2) * exp(I*6.2832*x) * exp(I*c*x^2).  Since exp(I*6.2832*x) is near 1 when x is an integer that's not too large,
we get a rather similar shape at intervals of 31416. 

Plot some partial sums of this series...

> S[0]:= 0;
   for i from 1 to 2*10^5 do S[i]:= evalf(S[i-1]+sin(.1e-3*i^2)) od:
   plots[pointplot]([seq([i,S[i]],i=1..2*10^5)],symbol=point);

Not at all what I might have expected to see.  Any explanation?

This is known as the "casus irreducibilis": for a cubic with rational coefficients and three real irrational roots, it is impossible to express the roots using only real radicals.  See e.g. en.wikipedia.org/wiki/Casus_irreducibilis.  On the other hand, you can express the roots using the functions sin and arctan.  For example:

> simplify(evalc([solve(x^3-5*x+1)]));

[2/3*5^(1/2)*3^(1/2)*sin(1/3*arctan(1/9*3^(1/2)*473^(1/2))+1/6*Pi), -1/3*5^(1/2)*(3^(1/2)*sin(1/3*arctan(1/9*3^(1/2)*473^(1/2))+1/6*Pi)+3*sin(-1/3*arctan(1/9*3^(1/2)*473^(1/2))+1/3*Pi)), -1/3*5^(1/2)*(3^(1/2)*sin(1/3*arctan(1/9*3^(1/2)*473^(1/2))+1/6*Pi)-3*sin(-1/3*arctan(1/9*3^(1/2)*473^(1/2))+1/3*Pi))]

Unfortunately the Maple type function is poorly named, it really should be something like "function_call", i.e. an expression of the form a(args) where a is a name and args is an expression sequence.  In your case I think you want y to be of type procedure.

The basic problem with what you're doing here is that the formal parameters x and z of E have no connection to the local variables x and z of the procedure elasticity.

I think you want something like this:

> elasticity:= proc(y:: procedure)
    (x,z) -> D(y)(x)*x/z;
  end proc;

  E := (x,z) -> elasticity(y)(x,z);