Why not just use @@ instead of ^ in the first place? 

MyOp:= ((A+B)@@2)(x); #your 1st question

subs(A= f, B= g, MyOp); #your 2nd question

Yes, the function is NumberTheory:-PrimeCounting. If your Maple is too old to have that, identical functionality is provided by numtheory[pi].

I suppose that your F is a PDF. If it is reasonably smooth^[1], then this can be done by Statistics:-Sample. See the third example on its help page, which deals with a custom distribution. I suspect that the efficiency of this example could've been improved by including the basepoints and range options. I'll look into that later.

^[1]Twice differentiable with finitely many inflection points on its support.

You are severely misreading the Wikipedia page that you linked to, which states clearly, in very prominent locations on the page, that A = L.L^* is the standard Cholesky factorization^[1] and that A = L.D.L^* is "a closely related variant."^[2].

Regarding "the output described in books": I think that you'd be hard pressed to find a more standard or respected textbook on the subject than Golub and Van Loan's Matrix Computations^[3], where is found

• Theorem 5.2-3: Cholesky Decomposition. If A ... is symmetric positive definite, then there exists a lower triangular G ... with positive diagonal entries such that A = G.G^T.

A pseudocode algorithm for computing G is given as Algorithm 5.2-1. The L.D.L^T decomposition is given as Theorem 5.1-2, but no name, Cholesky or otherwise, is attributed to it.

^[1]This is in the first sentence of the second section.

^[2]This is in the first sentence of the third section.

^[3]Gene H Golub and Charles F Van Loan, Matrix Computations, 1983, Johns Hopkins University Press, ISBN: 0-8018-3010-9

Here is another way to do it. I think that it's more general.

simplify((lhs-rhs)(eqq)) assuming positive;^[1]
`if`(%::`*`, select(depends, %, tau), %) = 0;^[2]
collect(%, indets(%, {specfunc(diff), typefunc(q[posint])}));^[3] 

^[1]simplify operates independently on the two sides of an equation. In particular, it doesn't divide out common factors. (lhs-rhs)(eqq) subtracts the right side of the equation from the left side, forming an expression rather than an equation. The assuming positive allows the square roots to be simplified.

^[2]If the resulting expression is a product, select the factors containing tau; the rest are discarded.

^[3]collect with respect to any derivatives and any functions of q. Note that these do not need to be explicitly listed.
 

Consider the univariate case: limit(f(x), x= a). There are only two paths along which x can travel to get to a, and we call those "from the left" and "from the right" (and those limits may be different). But in the multivariate case, there are an infinitude of paths along which (x,y) can travel to get to (a,b), and there may be many limits.  To find those, you need to parameterize the path. For example,

X:= a+t; Y:= b-t; limit(f(X,Y), t= 0, right);

Since tanh(p) ->  1 as p -> infinity, just replace tanh(p) by 1 in HeunG and you get Float(infinity) + Float(infinity)*I. Now, I'll be the first to admit that this technique for finding a limit is not 100% reliable. But it can be backed up with numeric and graphical evidence.

The divergence of the imaginary part is much slower than the divergence of the real part.

From what you've said in your two most-recent replies, I'm now convinced that this epsilon nonsense is the source of your troubles. Let's summarize what you're doing:

1. You make up a somewhat arbitrary small positive number epsilon. You've used 1e-10. You have to admit that there's no scientific basis for this particular value.
2. You compute tan(Pi/2+epsilon), thereby obtaining an arbitrary negative number of very large magnitude (approximately -1/epsilon).
3. After several pages of calculations, you divide some terms by this number to produce yet another arbitrary small factor of approximately -epsilon.
4. You claim to do this because the terms being divided are "neglible".
5. You claim that when you do this in a "similar model", "I'm getting the results I want."
6. Since a larger value of Digits (15) is giving you results that you don't "want", you ignore that truth and settle for the false security of using a smaller value of Digits (10). 

When summarized like that, it seems almost crazy. The terms are neglible because you're dividing them by a large magnitude; that doesn't mean that they were inherently neglible. If for some other valid reason they are in fact inherently neglible, why did you expend the effort to type them into the formula? You should simply omit them from the formula and include a comment such as "I am omiting the term ___ because it is neglible because of (some scientific reason)."

Are you under some political pressure to include those terms and you're trying to obfuscate the fact that they're artificially being removed?

While it's true that you can't divide by tan(Pi/2) in a numerical computation, you can multiply by cot(Pi/2) = 0. Isn't that equivalent?

A scientist's goal ought to be getting correct results, not the results that they want.

The differential order of the equation is 2, so you need 2 initial or boundary conditions, whereas you've supplied 1. The error message indicates that it's expecting you to supply a value for D(u)(0), i.e., u'(0). So, you can add something like D(u)(0) = 1 to SYS2.

It would help if you show your code. When you say "the numerical solver", I assume that you mean LinearAlgebra:-Eigenvectors. I'll call your constructed matrix M. As you no doubt realize, M is also symmetric. The solver can produce more accurate results if you tell it this. So, instead of doing

LinearAlgebra:-Eigenvectors(M),

do

LinearAlgebra:-Eigenvectors(Matrix(M, shape= symmetric, datatype= float[8])).

Without having your actual matrix, I can't be sure if this will help. So, if it doesn't, please post your actual matrix.

You can't have a derivative satisfy some nice, smooth, normal condition at all points except for a few specific points at which you specify its value. It makes no sense mathematically. You'd might as well say that you can have more than one boundary value per derivative.

The Answer by user sand15 shows you how to correct the syntax so that dsolve won't reject your input, but its output is not a solution to your posed problem, nor can it ever be. On the other hand, if you change the problem to include an epsilon, as in his Reply to Kitonum, then it becomes possible.

The computational difference that gc() makes is the deletion and reinitialization of the remember tables of all procedures declared with option system. I suspect that the additional anomaly that you're having with kernelopts(memusage) is because its usage is coincidentally provoking a call to gc().

So, it would be entirely appropriate for you to call gc() before simplify. The only reason that use of gc is discouraged is because of its potential impact on efficiency (i.e., each extra gc() takes some time).

Your second argument to int should be of type list(name= range), such as [px= -infinity..infinity, py= -infinity..infinity, pz= -infinity..infinity], which could be abbreviated to [px,py,pz]=~ -infinity..infinity.

I'm not claiming that this will give you the result that you ultimately want. I'm just saying that you need to do this to overcome the syntax error that you currently face.

This is a univariate problem with finite bounds and no other constraints. For such a problem, the algorithm "Quadratic Interpolation"---which doesn't use an initial point---may be used; and indeed it is used in this case. You can see which algorithm is used and get a bunch of other information by setting infolevel[Optimization]:= 5. The algorithms are described on the help page ?Optimization,General,Methods. If change the method, then an initial point can be used:

Optimization:-NLPSolve(sin(x)/x, x= 1..30, initialpoint= {x = 17}, method= branchandbound);

It seems that pdsolve is having trouble with your boundary condition at infinity. So, I changed infinity to a variable, Infinity, before passing the input system to pdsolve, and I applied limit(..., Infinity= infinity) to the output returned by pdsolve after it was returned. This seemed to work (no undefined in the output), but I don't have time right now to check the accuracy of those results. Perhaps you can do that and report back here?