@rzarouf You're welcome. You found a Maple bug, and I have submitted an SCR (bug report), so it is fixed in future versions.

@Carl Love The method is fine, but M^%H.M is not symmetric, so using shape=symmetric forces the offdiagonals equal. So you should omit that.

@Scot Gould Although @Rouben Rostamian has clarified that the any Robin boundaries can be solved, it seems that Maple finds sometimes the right answer consistent for the initial condition given and sometimes not. Your first solution was already compatible for the right-hand boundary (at 50 points it is pretty close on almost all the right half of the domain). So I knew to play around with the left-hand one, and used an analogy with the solutions that I was most familar with to figure out that sign change. I usually solve cases that evolve toward steady state, so I added my bias there. But maybe Maple is also consistently solving some classes of boundary conditions; it would help if we know more clearly what assumptions it is making. In this respect, @nm's suggestion to put in the initial solution in afterwards might be counterproductive; maybe Maple is not returning a solution because it knows it can't solve the case compatible with that initial conditions. I suspect very specific assumptions on signs of parameters at the pdsolve stage might be helpful.

I don't think the explanation on the Wikipedia page is that helpful; one needs to think of the relative signs to achieve the right physical outcome. I deliberately kept the fluxes the same sign to get to a steady-state, then adjusted the other signs. The general steady state solution is a straight line. If the slope of that is positive, the flux is negative and there will be the same heat flux coming in the right-hand side as going out the left-hand side at steady state. So if the initial fluxes are not equal but of the correct sign, then the other parts of the Robin boundary condition have to allow the slopes to evolve to be the same (assuming you want steady-state of course). In your case the fluxes evolve to be zero.

I didn't really want to use all that opcrap to extract the pieces and did try to do it without a loop first, but I also didn't want to spend to much time optimizing it. Make a procedure that determines whether it is an eigenvalue solution without an analytical solution for the eigenvalues might be worthwhile, but in that case I would probably solve it "by hand". Actually with Maple that is not too difficult. If you call pdsolve with HINT=f(t)*g(x), then extracting the ode for t and x and relating them to the eigenvalue is straightfoward. The second-order ode in x can be solved for {g(x),lambda} and if there is an analytical formula for lambda, Maple will find it; otherwise you have to find the eigenvalues numerically.

I should say, I almost never solve the heat/diffusion equation by this eigenvalue method, because the error functions and x/sqrt(k*t) usually found don't easily come out of it. Instead I find the solution by Laplace transforms. (To get the general solution containing x/sqrt(k*t) you need to use SimilaritySolutions.)

@Rouben Rostamian  Thanks for the clarification. (Yes, I generally expect solutions to evolve toward a steady state.)

@Teep "Do you think that it is possible to extract a general expression for the roots of f (or g) in terms of a and b?" No, if Maple could find that it would not have returned the RootOf, which is the "next best thing" when it can't find an explicit solution. Sometimes it is possible to manually rearrange things to help Maple out, but I think with a non-integer b in the exponent, it is likely an intractable problem. Even for the integer b case, polynomials of high degree do not usually have explicit symbolic solutions.

@Saalehorizontale As you expect, a RootOf(..., 0.35..0.36) or RootOf(...,index=2) corresponds to a single root. Unfortunately, as you already saw, RootOf without an index or range or label sometimes means any/all of the roots and sometimes means the first one. Solving a polynomial gives a RootOf that says all possible roots work, and can be found by allvalues, but evalf just produces one. Other commands, particularly those where a RootOf stands for an algebraic number in a field mean also the first one, e.g. evala(Minpoly(..)).

I didn't check the timings, but it's not unexpected that figuring out the actual numerical values of the roots is slow. On the other hand solving for polynomial roots is optimized compared to roots of more general equation, e.g., fsolve(polynomial in x,x,complex) will quickly return all roots.

@dharr I edited the above to include floor(root(evalf(n), 3)). In floor(root(n, 3)), root(n,3) just gives n^(1/3) and floor does all the work.

@GFY You can solve it by narrowing the ranges, with guidance of a plot. You could modify the equations so a__3=0,c__3=0 is not a solution by factoring and dividing through by power(s) of a__3 and c__3.

question11.mw

@GFY That was unexpected for me. Solve is really for symbolic solutions, and you have floats in your equations. I always formulate equations symbolically, and try solve. Then one can substitute values for the parameters, e.g., params:= {a=3.1, b=2.1 etc}, then eval(eqn, params);

@Aixleft math Typically I specify as many variables as equations, but if there is a solution with more, Maple can often find it.

@mz6687 You want to make assumptions on all variables that you want to be real:
simplify(conjugate(q)) assuming a::real, t::real etc. (or make all assumptions at the beginning of your worksheet)

evalc automatically assumes all variables are real, so another approach would be to separate lambda1 into real and imaginary parts and then use evalc, which is probably the approach I would take.

@nm Interesting. Note A.5 gives 5*A but 5.A gives "Error, missing operator or `;`", and A.B gives A.B (the matrix multiplication symbol).

@mz6687 OK, but without further information the question remains as to why they should both give zero for the same set of parameters - that seems to be your expectation and question.
 PS. I edited my above comment since I forgot evalc will also assume lambda1 real. Probably using conjugate with assumptions on variables is the way to go. 

Just a couple of comments, based on what I think you want to do. You manually do a complex conjugate by substituting I=-I, lambda1 = conjugate(lambda1), etc, but simplify(qb-conjugate(q)) does not give 0. Edit: removed incorrect solution to this.

You have two equations

eq1 := -2*q*diff(qb,t) - 2*(diff(q,t))*qb + diff(v,x) =0;
eq2 := 2*I*diff(q, x,t) + I*q*v + I*v*q + 2*diff(q,t)=0;

and substituting in some parameters you find eq2 gives 0=0 but eq1 does not. But they look very different to me. Are they supposed to be complex conjugates of each other?

You have many subscripted Energy variables so I can't really understand what you want to do, and some change of variables is mentioned but without being specific. It seems strange to differentiate and then integrate almost the same thing. It would be a lot easier to understand if you have variables instead of numbers (see my example below).

Perhaps you could just solve a differential equation. If you just want plots, something like the following implicitly dies a change of variables from t to E. Or you could use PDEtools:-dchange to get Maple to do the change of variables on the original ode. Perhaps the following does something that can be adapted to your situation.

Download DAE.mw