@MaPal93 I can imagine some pathological cases, but in general I would expect that the results should be the same. If I use limit(X[i], gamma = infinity) in generic_n-form_of_infinity_limit2_MaPal.mw I don't get the expression for X[i] in  n_from_2_to_6_-_not_coincide.mw so perhaps something is incorrectly entered, or perhaps I am not understanding what you are doing.

@MaPal93 I rewrote with as few summations as possible and was able to get the following simple result:

generic_n-form_of_infinity_limit2.mw

The _a indices arise when it would otherwise lead to a summation over i inside another summation over i, which is confusing since they are only dummy indices and not related. So Maple makes sure they have different names.

@MaPal93 I don't have any specific suggestions other than try to figure out the patterns.

Download patterns.mw

@MaPal93 Actually, I realized that if you use collect(W, gamma) (first method) you see W = (stuff1)*gamma+(stuff2), so the limit will be +infinity if stuff1 > 0 and -infinity if stuff1 < 0. So it seems correct that is infinite. (But now you have new X[i] etc.)

In general if you see something like signum(n*rho - rho + 1) you can just make some assumption like n*rho - rho + 1>0 to make it go away, assuming that such an assumption is valid.

@salim-barzani So if the input is 1/G(xi)*($$$$$$$$$$$$+diff(G(xi),xi)*(##########) what is the output you want?

I don't understand what you are trying to do, but this might be useful.

collect.mw

I don't think there is any simple way. Since you want several matrices, GenerateMatrix won't work directly. You can use dcoeffs to find the pieces and then assemble them into the matrices.

@janhardo The DLMF (and the corresponding hardback book version) is produced by NIST (National Institute of Standards and Technology) in the U,S., and is the successor to Abramowitz and Stegun's, "Handbook of Mathematical Functions"; it is considered an authoritative source, perhaps "the" authoritative source.

These differential equations are standard forms in the sense that they are often quoted, e.g., one speaks of "Bessel's equation"  and there is a common way to write it, though of course there is not universal consistency in these things. 

But in terms of classifications in a text as I think you want, there is probably not much of a standard. I personally like Zwillinger's "Handbook of Differential Equations". Of course these classifications exist because of different solving methods, but since the objective is often a special function, they are related to the other standard forms.

Maple' odeadvisor(verg) gives: [[_2nd_order, _missing_x], [_2nd_order, _reducible, _mu_x_y1]]

Since many of the second-order differential equations and their "special" functions are related through transformations to the Sturm-Liouville eqation, perhaps it is the standard form to rule them all :-).

@MichaelVio Here's an example.

Download dchange.mw

@salim-barzani I'm not sure I fully understand the order you want to do things or what you want to store, but perhaps this:

Download find_all_case_solution_of_ode_.mw

For Vectors, use a:=<2,1>, not <<2,1>>, which is a 2x1 Matrix. Then you can use "." for the dot product.

@Carl Love I think shape = symmetrix or shape = hermitian change the algorithm for float matrices to produce only real eigenvalues and eigenvectors but have no effect in the symbolic case.

@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.)