As you read through this explanation, if you get to any point that you don't understand, please ask about it. You need to understand each step before going to the next step.

First you should strive to understand s(x). Using some plots if you need to, convince yourself that:

1. s is periodic with period 1,
2. on its fundamental period 0 <= x <= 1, s is equivalent to piecewise(x <= 1/2, x, x <= 1,1-x).

Let me know if you need help with the plot command.

It follows that s'(x) is also periodic with period 1, and is equivalent on its fundamental period to 

piecewise(x=0, undefined, x < 1/2, 1, x = 1/2, undefined, x < 1, -1, x = 1, undefined). 

Note that s'(x) is always an integer when it's defined. We'll use that fact later. Also, s'(x) is defined for all real x except when x = p/2 for some integer p.

Now define S_M(x) = s(M*x)/M, for any positive integer M. It follows from basic facts about periodic functions (no calculus required) that S_M(x) has fundamental period 1/M. From the chain rule, we have that S_M'(x) =  s'(M*x) wherever it's defined. In particular, it's either 1 or -1 wherever it's defined. And it's defined for all real x except when x=p/(2*M) for some integer p.

Proposition 1: If you have two periodic functions f(x) and g(x) with fundamental periods P and Q and P/Q is an integer, then the period of h(x) = f(x) + g(x) is P. 

[Edit: The red characters have been deleted from the next line.]

Now define b_N(x) = sum(s(2^n*x)/2^n, n= 1..N). Its period is 1/2^N (by that proposition) and it's differentiable for all real x except when x = p/2^(N+1) for some integer p. When it's defined, b_N'(x) is a sum of exactly N values, each being either 1 or -1. So, it's an integer between -N and N, inclusive.

Your original problem was b₄'(9.05), which equals b₄'(.05). The critical points are p/32 for any integer p. The value is S₂'(.05) + S₄'(.05) + S₈'(.05) + S₁₆'(.05) = 1 + 1 + 1 + (-1) = 3 because 1/32 < .05 < 2/32.

Note that we can easily compute b_N'(x) for any positive integer N and any real x by looking at the first N digits after the decimal point in the binary (base-2) expansion of frac(x). You just add 1 for every 0 and subtract 1 for every 1.

It can be created easily like this:

n:= 4:
M:= Matrix(n, scan= band[1,0], [[1$(n-1)], [$1..n]]);

So, the band[1,0] says that there is 1 subdiagonal and 0 superdiagonals, and the initializer is given by those diagonals, in that order.

To optimize the storage, use

M:= Matrix(
   n, 
   shape= triangular[lower], storage= band[1,0],
   scan= band[1,0], [[1$(n-1)], [$1..n]]
);

So, when you use scan= band[...], then you can make your initializer logically follow your diagonal structure, rather than using some complicated indexing scheme that goes by rows and columns.

Yes, it is easy to build a matrix by its submatrices. I usually use the angle-bracket constructors for that. If you give an example, I'll create it.
 

I implemented the L.D.L^T algorithm for hardware-float matrices for you. It's not as fast as NAG code, but it's not bad. I think that you'll find it useful. Let me know.

Would someone who knows more about compiled Maple than I do suggest some improvements for the timing of this?
 

Download LDLT.mw

The seq command does preserve the order. However, you took the output of seq and made it into set by using {...}. Sets do not preserve order, which is consistent with their mathematical definition. If you simply remove those {...}, then you'll get what you expect.

You can assign a default value to a parameter, and thereby make it optional, by directly assigning to it in the proc line:

sim:= proc(x, n, o:= exactly) 

You can do this and also put a type restriction on the parameter like this:

sim:= proc(x, n, o::identical(exactly,minimum,maximum):= exactly)

There's a huge number of different parameter types, and those help pages that you referenced could use a lot more examples.

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.