This is in reference to your original system in which you say "PDESolve function gives a very strange solution".  pdetest verifies the solution: pdetest(sol, [PDE, bc1, bc2, bc3]) gives [0, 0, 0, 0]; also easily seen by hand. Since you are setting the Laplacian equal to zero, you can omit the (cosh(eta) - cos(xi))^2 factor and just use diff(u(xi, eta), xi, xi) + diff(u(xi, eta), eta, eta) = 0. So the problem is the same as for a rectangular region in which there is no variation in the xi direction. In other words it reduces the 1-D diffusion over a finite range, which does just have the linear solution.

If you upload your updated worksheet, I'll take a look.

p.s. you can use PDE := VectorCalculus:-Laplacian(u(xi, eta), 'bipolar'[xi, eta]);

@GFY You now mention a stability analysis, which IMO is more easily done via the Jacobian than Liapunov. In the attached I watch the eigenvalues move as a function of U, and evidence is seen of a Hopf bifurcation. (The complexplots do not seem that accurate and could probably be optimized.) Not sure if this is useful to you; you seem more interested in other equilibria, which I addressed in my earlier reply.

As I said before I don't see how the Bifurcation command can be useful here, but maybe someone else can think of a way.

Download bifurcation.mw

The bifurcation command is only for iterative maps, not systems of odes. If I understand correctly, you want to know how many equilibrium points (is this what you mean by solutions?), stable or not, that exist for various values of the parameter U. Is this correct?

To find this I would convert your system to a system of six first-order odes dQ_i/dt = f_i(Q_j), and then find the number of solutions to the polynomial system f_i(Q_j) = 0, i=1..6. I don't think you can reliably do that with fsolve and floating point coefficients as you have. So I would change the floating point numbers to symbolic parameters with the smallest number of independent parameters. Even then the system may be too complicated to work with, but if you can simplify it in this way I'll take a closer look.

@sand15 Excellent! Indeed, there seems to be no ternary plot yet.

@Carl Love Sorry; I realized that afterward. As you say, the is() returning false in the original question remains puzzling, though I guess in some sense it is false for n = 2.5. Since the product in term2 doesn't evaluate to any simpler expression, it's hard to see why it doesn't always return FAIL.

@AHSAN As I said, "So I expect u(y) to go to 1 at y=1 and to have a slope of zero at y=0, both of which seem to be true from your plots." Perhaps I don't understand what you mean by start at zero - for me that means that u = 0 at y=0. But depending on the curve, u is in the range 0.45..0.5 at y=0, which is not zero and not 1.

@AHSAN For u(y), which is D(psi)(y) you have

D(psi)(1 + x^2/2) = 1, D(D(psi))(0) = 0

So I expect u(y) to go to 1 at y=1 and to have a slope of zero at y=0, both of which seem to be true from your plots.

@AHSAN You say theta(y) and phi(y) should start from 1, but you have

D(theta)(1 + x^2/2) = 1, theta(0) = 0,
D(phi)(1 + x^2/2) = 1, phi(0) = 0

so I expect theta and phi to start from zero and have a slope of 1 at y = 1 (since x=0).

I don't use Safari, so can't comment about browser issues. But if I go into my account information, I can change my email address and keep the same account. So I suggest you ask technical support to change your account to another email address. Then logging on with that email address and using the forgotten password mechanism should solve that problem.

@Math-dashti

 true.mw

@mehdi jafari OK, see @mmcdara's answer for why this is hard to do. (I was confused by your terminology since I do not consider finding where FF(x,P1) = 0 to be "plotting the function FF".)

@Math-dashti Like this?

Download solve.mw

@mmcdara The worksheet opens with the comma's but executing it changes them to ".", but the plot error is the same owing to complex values. 

@Math-dashti I don't understand exactly what you want - please specify in the notation of your worksheet exactly what you want. I think S__1 and S__2 are some functions of p1,p2,p3,k? If so I don't know what  functions those are and what you want to do with them?

Or if you want to reproduce the paper results, please give a worksheet with that notation.

@Andiguys As I said, the view option is the syntax to magnify it. but it isn't helpful in this case.