According to the OEIS entry, there is a "closed form", if you count 2F1 hypergeometric functions:

p(n) = hypergeom([2-n, n+1], [2], -1)

Just to clarify: what that "real" means is not that f is real, it is that x is real.  So for example

> limit(exp(-1/x^2), x=0);

0

(even though the limit as x -> 0 in the complex plane would not exist).

The problem (which still exists in Maple 12) seems to have to do with the nested RootOfs.  For a simpler example:

> Q:= RootOf((p-RootOf(_Z^3-1))*_Z^3-1);

> eval(Q, p = 2.0);

Error, (in evala/preproc3) floats not handled yet

Technically, this may be a weakness rather than a bug.  I think there are some difficult mathematical issues connected with the interaction between algebraic numbers and floats.

A work-around is to either use exact rationals instead of floats (e.g 22/10 instead of 2.2), or to substitute the parameter values before trying solve.  The second is probably better.  Thus:

> S:=solve(eval({op(F),J[2,2]},pars),{x,y,A});
    remove(has,{S},I);

{{A = -.9791368510,x=-0.7190485833,y=0.1193582308},{A=0.1380471481,x=0.01107628534,y=1.256721040}}

The problem (which still exists in Maple 12) seems to have to do with the nested RootOfs.  For a simpler example:

> Q:= RootOf((p-RootOf(_Z^3-1))*_Z^3-1);

> eval(Q, p = 2.0);

Error, (in evala/preproc3) floats not handled yet

Technically, this may be a weakness rather than a bug.  I think there are some difficult mathematical issues connected with the interaction between algebraic numbers and floats.

A work-around is to either use exact rationals instead of floats (e.g 22/10 instead of 2.2), or to substitute the parameter values before trying solve.  The second is probably better.  Thus:

> S:=solve(eval({op(F),J[2,2]},pars),{x,y,A});
    remove(has,{S},I);

{{A = -.9791368510,x=-0.7190485833,y=0.1193582308},{A=0.1380471481,x=0.01107628534,y=1.256721040}}

I didn't get that bug in Maple 9.5.  Perhaps you could upload your worksheet.  I did get something else interesting: two spurious solutions with y=0, which clearly don't satisfy J[2,2] = 0.  Disregarding those, I get three bifurcations:

[{a = -.9628740975, y = .1190550790, x = -.8162589472}, {a = 1.249060056, x = .1592138727, y = .1190550789}, {y = .1190550789, a = .9588448130, x = .8747161127}]

Here's my worksheet.

Download 4541_bifurcation.mws
View file details

I didn't get that bug in Maple 9.5.  Perhaps you could upload your worksheet.  I did get something else interesting: two spurious solutions with y=0, which clearly don't satisfy J[2,2] = 0.  Disregarding those, I get three bifurcations:

[{a = -.9628740975, y = .1190550790, x = -.8162589472}, {a = 1.249060056, x = .1592138727, y = .1190550789}, {y = .1190550789, a = .9588448130, x = .8747161127}]

Here's my worksheet.

Download 4541_bifurcation.mws
View file details

1)  It looks like a bug.  Can you give us the exact commands you are using?
I don't know what are F2, E2 etc.

2) Yes, it should be possible.

1)  It looks like a bug.  Can you give us the exact commands you are using?
I don't know what are F2, E2 etc.

2) Yes, it should be possible.

1) I did not use the initial conditions, because in my understanding of bifurcation theory initial conditions are not important.  Bifurcations occur when equilibria (or sometimes other invariant sets) change their character. 

2) The bifurcation occurs at approximately e = 85.43224375.  In order to see what happens there, I looked at one value of e on each side of this point. 

1) I did not use the initial conditions, because in my understanding of bifurcation theory initial conditions are not important.  Bifurcations occur when equilibria (or sometimes other invariant sets) change their character. 

2) The bifurcation occurs at approximately e = 85.43224375.  In order to see what happens there, I looked at one value of e on each side of this point. 

Maple 9.5 needs sets rather than lists in solve.  So it should be

 S:= solve({op(F), J[2,2]}, {x,y,e});

Also, Maple 9.5 does not allow the gridrefine option in implicitplot.

Using J[1,1] will not work, because F[1] and J[1,1] do not depend on y or e.

_Z is used as the variable in RootOf notation.  That is, RootOf(f(_Z)) stands for any value _Z such that f(_Z) = 0.

Maple 9.5 needs sets rather than lists in solve.  So it should be

 S:= solve({op(F), J[2,2]}, {x,y,e});

Also, Maple 9.5 does not allow the gridrefine option in implicitplot.

Using J[1,1] will not work, because F[1] and J[1,1] do not depend on y or e.

_Z is used as the variable in RootOf notation.  That is, RootOf(f(_Z)) stands for any value _Z such that f(_Z) = 0.

The right side of the system:

> F := [a*b*((x^2+c^2)/(x^2+1))-d*x, a*e*((y^3+f^3)/(x^3+1)*(g*x+1))-y];

The Jacobian:

> J := VectorCalculus[Jacobian](F, [x,y]);

J := RTABLE(45781648,MATRIX([[2*a*b*x/(x^2+1)-2*a*b*(x^2+c^2)/(x^2+1)^2*x-d, 0], [-3*a*e*(y^3+f^3)/(x^3+1)^2*(g*x+1)*x^2+a*e*(y^3+f^3)/(x^3+1)*g, 3*a*e*y^2/(x^3+1)*(g*x+1)-1]]),Matrix)

A bifurcation can occur when the real part of an eigenvalue of the Jacobian at an equilibrium point crosses 0.  In this case the Jacobian is lower triangular, so the eigenvalues are the diagonal elements (and are real).

> S:= solve([op(F), J[2,2]], [x,y,e]);

 S:=[[x = RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z), y = RootOf(2*_Z^3-1,label = _L1)*f, e = 2/3*(d+a*b*RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z)^2+a*b*c^2-d*RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z))/d/a/f^2*RootOf(2*_Z^3-1,label = _L1)/(g*RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z)+1)]]

 Try it for some values of the other parameters:

 > pars:= { a=7, b = 0.2, c=0.2,d=0.1,f=1,g=0.1};

> evalf([allvalues(eval(S,pars))]);

[[[x = 13.93110365, y = .7937005260, e = 85.43224375]], [[x = .3444817452e-1+.1975123882*I, y = .7937005260, e = .7499096375e-1-.2003550007e-2*I]], [[x = .3444817452e-1-.1975123882*I, y = .7937005260, e = .7499096375e-1+.2003550007e-2*I]], [[x = 13.93110365, y = -.3968502630+.6873648188*I, e = -42.71612188+73.98649342*I]], [[x = .3444817452e-1+.1975123882*I, y = -.3968502630+.6873648188*I, e = -.3576035666e-1+.6594585468e-1*I]], [[x = .3444817452e-1-.1975123882*I, y = -.3968502630+.6873648188*I, e = -.3923060707e-1+.6394230467e-1*I]], [[x = 13.93110365, y = -.3968502630-.6873648188*I, e = -42.71612188-73.98649342*I]], [[x = .3444817452e-1+.1975123882*I, y = -.3968502630-.6873648188*I, e = -.3923060707e-1-.6394230467e-1*I]], [[x = .3444817452e-1-.1975123882*I, y = -.3968502630-.6873648188*I, e = -.3576035666e-1-.6594585468e-1*I]]]

We're only interested in real values of x,y and e, so the bifurcation should be at e= 85.43224375.

For e = 85.4:

> F854:= eval(F,{op(pars),e=85.4});
   with(plots): with(DEtools):
   display([implicitplot(F854,x=13.92 .. 13.95,y=0.7 ..0.9,gridrefine=2,colour=[green,blue]),
   DEplot(zip(`=`,[D(x)(t),D(y)(t)],eval(F854,{x=x(t),y=y(t)})),[x(t),y(t)],t=1..20,x=13.92 .. 13.95, y=0.7 .. 0.9)]);

Here we see two equilibrium points, one a saddle and the other a stable node.

For e = 85.5:

> F855:= eval(F,{op(pars),e=85.5});
  display([implicitplot(F855,x=13.92 .. 13.95,y=0.7 ..0.9,gridrefine=2,colour=[green,blue]),
  DEplot(zip(`=`,[D(x)(t),D(y)(t)],eval(F855,{x=x(t),y=y(t)})),[x(t),y(t)],t=1..20,x=13.92 .. 13.95, y=0.7 .. 0.9)]);

The two equilibrium points have collided and disappeared.

The right side of the system:

> F := [a*b*((x^2+c^2)/(x^2+1))-d*x, a*e*((y^3+f^3)/(x^3+1)*(g*x+1))-y];

The Jacobian:

> J := VectorCalculus[Jacobian](F, [x,y]);

J := RTABLE(45781648,MATRIX([[2*a*b*x/(x^2+1)-2*a*b*(x^2+c^2)/(x^2+1)^2*x-d, 0], [-3*a*e*(y^3+f^3)/(x^3+1)^2*(g*x+1)*x^2+a*e*(y^3+f^3)/(x^3+1)*g, 3*a*e*y^2/(x^3+1)*(g*x+1)-1]]),Matrix)

A bifurcation can occur when the real part of an eigenvalue of the Jacobian at an equilibrium point crosses 0.  In this case the Jacobian is lower triangular, so the eigenvalues are the diagonal elements (and are real).

> S:= solve([op(F), J[2,2]], [x,y,e]);

 S:=[[x = RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z), y = RootOf(2*_Z^3-1,label = _L1)*f, e = 2/3*(d+a*b*RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z)^2+a*b*c^2-d*RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z))/d/a/f^2*RootOf(2*_Z^3-1,label = _L1)/(g*RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z)+1)]]

 Try it for some values of the other parameters:

 > pars:= { a=7, b = 0.2, c=0.2,d=0.1,f=1,g=0.1};

> evalf([allvalues(eval(S,pars))]);

[[[x = 13.93110365, y = .7937005260, e = 85.43224375]], [[x = .3444817452e-1+.1975123882*I, y = .7937005260, e = .7499096375e-1-.2003550007e-2*I]], [[x = .3444817452e-1-.1975123882*I, y = .7937005260, e = .7499096375e-1+.2003550007e-2*I]], [[x = 13.93110365, y = -.3968502630+.6873648188*I, e = -42.71612188+73.98649342*I]], [[x = .3444817452e-1+.1975123882*I, y = -.3968502630+.6873648188*I, e = -.3576035666e-1+.6594585468e-1*I]], [[x = .3444817452e-1-.1975123882*I, y = -.3968502630+.6873648188*I, e = -.3923060707e-1+.6394230467e-1*I]], [[x = 13.93110365, y = -.3968502630-.6873648188*I, e = -42.71612188-73.98649342*I]], [[x = .3444817452e-1+.1975123882*I, y = -.3968502630-.6873648188*I, e = -.3923060707e-1-.6394230467e-1*I]], [[x = .3444817452e-1-.1975123882*I, y = -.3968502630-.6873648188*I, e = -.3576035666e-1-.6594585468e-1*I]]]

We're only interested in real values of x,y and e, so the bifurcation should be at e= 85.43224375.

For e = 85.4:

> F854:= eval(F,{op(pars),e=85.4});
   with(plots): with(DEtools):
   display([implicitplot(F854,x=13.92 .. 13.95,y=0.7 ..0.9,gridrefine=2,colour=[green,blue]),
   DEplot(zip(`=`,[D(x)(t),D(y)(t)],eval(F854,{x=x(t),y=y(t)})),[x(t),y(t)],t=1..20,x=13.92 .. 13.95, y=0.7 .. 0.9)]);

Here we see two equilibrium points, one a saddle and the other a stable node.

For e = 85.5:

> F855:= eval(F,{op(pars),e=85.5});
  display([implicitplot(F855,x=13.92 .. 13.95,y=0.7 ..0.9,gridrefine=2,colour=[green,blue]),
  DEplot(zip(`=`,[D(x)(t),D(y)(t)],eval(F855,{x=x(t),y=y(t)})),[x(t),y(t)],t=1..20,x=13.92 .. 13.95, y=0.7 .. 0.9)]);

The two equilibrium points have collided and disappeared.

showstat(`evalf/exp`) is not that enlightening.  For floats, at sufficiently low settings of Digits, Maple just passes the computation off to evalhf, which uses the floating-point hardware of the system.  To see how exp is implemented there, you'll have to ask Intel, not Maple.  If Digits is too high to use evalhf, you want to look at `evalf/exp/general`.