@Alfred_F The last line of Carl's code does this, but only for a given value of k. So

combine(sin(3*x)*cos(x)^3)

gives

1/8*sin(6*x)+3/8*sin(4*x)+3/8*sin(2*x)

As usual you can go the other way with expand. But few Maple commands will deal with the general case; diff is an exception. It is this issue that makes your example hard. For example

int(sin(3*x)*cos(x)^3/x, x=0..infinity);

immediately gives 7*Pi/16.

@salim-barzani For indefinite integration there is always an arbitrary constant (which Maple doesn't show). If the antiderivative is F(x), then indefinite untegration gives F(x) + c, or just F(x) in Maple. Definite integration from -infinity to x gives F(x) - F(-infinity) and the assumption is that F(-infinity) = 0 since we haven't accumulated any area yet at the start. Maple has Intat to describe integration at a point - you could use Intat(f(a),a=x) instead of indefinite integration.

The main Maplesoft site has said "site under maintenance" with possibly reduced functionality for about the same period (still true today), so looks like the revamp took Mapleprimes offline. The application center was also down for at least some of it.

@nm @Alfred_F Since odetest gives 0 for y(0)=a, it is probably doing something simple like below, though I didn't check its code.

dsolve.mw (won't display right now)

@Alfred_F Maple agrees with you.

expand(PDEtools:-dchange(y(t)=sqrt(u(t)),-ode/t,[u(t)]));
DEtools:-odeadvisor(%);

@acer Thanks for the vote and tips, especially for the limit. I knew there would be another answer for the labelling issue - I tried a few things with typeset and uneval quotes, before moving on.

I've become sold on the advantages of equations that are expressions implicitly equal to zero (I know @nm holds the opposite view). Then there aren't things to cancel, just things to divide by. Many things don't map across the `=` and I find normal and numer etc convenient in this context. It has the disadvantage that it is less readable, but I can live with that.

I solve this problem in Maple in a post here.

Incidently, evalf(limit(...)) finds the wrong value of the limit - seems to find the maximum point. I've submitted an SCR.

Download limit.mw

@MichalKvasnicka Here is a routine to find the matrix exponential by solving the ode system dX/dt = M.X. I haven't optimized it since it seems you might be wanting to port to another language. If you want it further optimized, I can do that, especially through the parameter mechanism, where you only find the procedure to solve the odes once, and then evaluate it for multiple different parameter sets.
[Edit: now calls dsolve only once.]

Download ExpByODEs.mw

@MichalKvasnicka So at the end you want to supply numerical values of the parameters and output a numerical matrix, is that correct?

In general the methods Maple chooses for symbolic solutions are different from those for numerical calculations, and putting numbers in at the end may give numerical errors that would not arise from methods optimized for numerical issues. This is especially the case for matrix exponentials. An instructive read is "Nineteen Dubious Ways to Compute the Exponential of a Matrix", doi: 10.1137/S00361445024180 or here. The summary is unless you have a normal matrix there are no general reliable methods, and special difficulties can arise if your matrix is non-diagonalizable (defective; missing eigenvectors) as is true for your case.

If you continue with your strategy and you use Maple at the end of the day, experiment with increasing Digits to make sure you have convergence.

Conceptually, the ode solution method is probably the easiest to implement (with a stiff method - in Maple use stiff=true).

But since the need for the matrix exponential come about from a way to solve some linear equations, can you solve them directly without needing the matrix exponential at all?

For more help, it would be more useful to know more specifically what you intend to do.

@acer  - thanks for your question that prompted the clarification - I had assumed it was for compact display, such as would be required in a description of a method or proof.

@MichalKvasnicka That's a good idea, but you have many variables, so I'm not sure if this meets your definition of compact.

In answer to your questions, you can just pass the whole Matrix:

optimize(expFT, 'tryhard')

@MichalKvasnicka In my experience, it is always better to remove dependent variables early. (Non-dimensionalization is one way to do this.) For example, compare 

CharacteristicPolynomial(F, x);
CharacteristicPolynomial(F2, x);

and you will see that 

- 

has been simplified to 1 in the x^2 coefficient. I chose to leave omega__2 in and not omega__1 because some other expressions were simpler. Of course you may lose some symmetry this way, but I would argue that extra variables are more likely to be a problem - also computationally this speeds things up.

You can try other substitutions, and I would be guided by nondimensionalization in selecting these.

In some commands you can use implicit=true, and get your answer in terms of RootOfs, which can then be given a symbolic label, but this is not available for MatrixExponential.

You can try algsubs, e.g. algsubs(P=above expression, F) but the type of expressions it can handle is limited, e.g., it won't do things with sqrts in.

subs will sometimes work, e.g. subs(f__d2*lambda__2*omega__2 + f__d1*lambda__1 + f__p*lambda__1 + f__p*lambda__2 = P, g) does work where g is the output of CharacteristicPolynomial(F, x); But this depends on how the expression is represented natively, which is not always the way it appears to be, so one needs to have a good knowledge of Maple to do this.

A major problem is picking out a subexpression amongst many. subsindets can help with this, but is more advanced Maple; I'm not sure of your expertise here. Some people like pattern matching, which I don't have experience in.

So the real answer to the question is to understand what Maple is doing in terms of the algorithm. If your matrix were diagonalizable one can formulate it in terms of exps of the eigenvalues, from which the awkward square roots arise. But we are short one eigenvector, so the exp is done through the Jordon form and the primary matrix function, involving derivatives evaluated at the eigenvalues. That seemed too complcated, so I didn't pursue it.

I'd say this sort of thing is in general very difficult, as you supposed.

Numerical solution shows, for several different a values, that the [1,1] point is OK as an initial condition; there doesn't seem to be any singularity or other issue with it.

Download internal_error_instead_of_no_solution.mw

@salim-barzani If I understand you correctly, you are asking where the transformation equations for xp, yp come from.

Download Y-pde-line2.mw

@Alfred_F In general, Maple is not going to give information that will verify an output (there are other computer systems for proofs in which every step is justified). Setting infolevel may give information, as @nm showed for the integral in the other thread. Here infolevel in the one-line limit that @Scot Gould showed gives no information. As he showed, you can also look into the actual code used, though this is a painful and often fruitless endeavor.

But for limits one can imagine that asymptotic series is likely to be a useful method. So one could do it partly by hand as below. Of course, then can then ask how Maple got the series, but you probably know the mechanics of this.

Download limit.mw