@salim-barzani Here's an attempt that circumvents the problem, but I would guess that your H is not actually a solution of the pde.

n1.mw

Edit - I really don't understand what is happening to the integrals there; in particular how the integration constants are dealt with in a consistent way.

@Alex0099 As far as I know, Maple does not have specific support for Green's functions. One just solves the equations with Dirac, as you did. You may or may not get a solution, depending on how difficult the problem is. Here are some examples in cartesian coordinates, which illustrate varying degrees of success.

Download Dirac.mw

Do you have a functional form of source term in mind for which the Green's function is just an intermediate step, or is the Green's function itself the objective? If you can bypass it that might be easiest.

@Alex0099 The error is because you used floating point 1.1 and 0.5 instead of 11/100 and 1/2. But since you want a symbolic solution, I would leave these as eta1 and eta2.

I'm not clear how much you want to help Maple along. In general pdsolve is not good with ICs/BCs and I would seek a general solution first and then apply the ICs/BCs afterward. It looks like it doesn't find a general solution with the source term. You could try more specific hints, but in the absence of that I would probably follow along with the procedure you gave in the earlier worksheet, but for the specify sinh(xi) etc BCs.

bc4 still needs some work. You want the integrated flux oround the two perimeters equal? In any case you should use D and not diff for the bcs.

PDE_upd_2.mw

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