@vv The OP's point is that the L.D.L^T decomposition can be computed without first computing the L.L^T decomposition, and that this can be done without using square roots. This is described in the Wikipedia page.

@Klausklabauter I'll see if I can implement the square-root-free L.D.L^T algorithm in Maple for you, guided by the Wikipedia article. Are you only doing the hardware-float case? Note that avoiding square roots is not necessarily such a great thing because it tends to be implemented in hardware these days, and thus is fast. Here's an example comparing sqrt on a million elements versus squaring those same elements, both in hardware floats:

M:= LinearAlgebra:-RandomMatrix(1000$2, generator= rand(0. .. 1.)):
M2:= copy(M):
CodeTools:-Usage(map[evalhf,inplace](sqrt, M2)):
memory used=7.63MiB, alloc change=7.63MiB, cpu time=125.00ms, real time=112.00ms, gc time=0ns
CodeTools:-Usage(map[evalhf,inplace](`^`, M, 2)):
memory used=7.63MiB, alloc change=7.63MiB, cpu time=141.00ms, real time=148.00ms, gc time=0ns

The point of my earlier messages was not that the LDL^T decomposition was somehow inferior to LL^T; my point was simply that it is not standard to call it "Cholesky".

@Mariusz Iwaniuk Vote up. Conversion to Heaviside is another way to simplify piecewise's conditions. If you convert the Heaviside expression back to piecewise, it'll be in the simplified form:

convert(h1, piecewise, t) assuming n::posint;

@Rouben Rostamian  Vote up. It is curious that the following does not work: 

inttrans:-fourier(piecewise(abs(t) < 2*Pi*n, cos(t), 0), t, w) assuming t::real, n::posint;

I think that piecewise often works better when the conditions are simplified, as you did, even if that means there are more conditons.

@Klausklabauter When your instructor told you not to use the "old method", I think that they are referring to the long-since deprecated Maple command linalg[cholesky]. I don't think that they are referring to any distinction between L.L^T and L.D.L^T.

@mmcdara That's what I thought. The pop-up is not saying that the input is incorrect. It's saying that there are two possible interpretations, and it's giving you the opportunity to select one of them. So, I stand by my statement that the OP's input is correct when entered as 2D Input, even in Maple 2015.2.

I also abhor the 2D Input, and this dialog box is one of the many reasons. I think that they managed to get rid of it for Maple 2018. There's some connection between that and the option function_assign that you see when the procedure is copied to a 1D Input field.

@Britzel Yes, that's what it's supposed to mean. If there's a case where that's not true, it's a bug. (Of course, returning undefined isn't considered as returning "a value" for the purposes of this question.)

@Kitonum You'd might as well use

signum(a - b)  assuming a < b;

because by your method, the expression is only being analyzed at this superficial level. How is the user expected to know which parts need to be encapsulated, and, if they did know that, why would they use Maple for this?

Is the third problem in your PDF simply to be ignored? It is completely different from the first two. Indeed, it's not even a numeric problem as no initial condition is given.

@Carl Love After reading the example shown by dharr, I realized that each of the eigenvectors corresponding to eigenvalue 0 has with itself precisely the type of symmetry that you want. These eigenvectors are always of the form < x, -x >.  Let < u, v > = < x, -x >. Any scalar multiple of an eigenvector is also an eigenvector for the same eigenvalue, so - < x, -x > = < -x, x > = < v, u >, corresponding to eigenvalue -0 = 0.

This is fairly easy to prove rigorously. (I'd consider giving it as an exercise to an honors undergraduate linear algebra class.) Let me know if you need help with that.

So, the symmetry anomaly that you're seeing is not related to the numerical anomaly. And I think that the numerical anomaly will resolve when you specify shape= symmetric.

@jefryyhalim Yes, I can reproduce the erratic behavior as you described. My first reaction is to increase Digits to 30. I start getting undefineds on line sol_L:=.... By the maxim I just stated, the results for all lower values of Digits are unreliable.

See how far you can get doing exact computation only. Convert all decimal numbers in the input into their equivalent fractions. Now, this may quickly get bogged down with expression swell. Or it may work. I recommend that you execute line by line. Remove all colon statement terminators so that you can observe the output of every line. As soon as you see a decimal anywhere in the output, correct the corresponding input to exact values.

I am suggesting this, at the moment, merely as a means to find the error (if any) in your logic. I'm not suggesting, at the moment, that the final production code will need to be done in this exact form.

@jefryyhalim Here's a maxim for all numeric computation (not just Maple): Never ignore the results provided by using a higher value of Digits when they differ strangely or significantly from those provided by a lower value unless an expert in numeric computation has looked over your problem and given you a specific computational/mathematical reason to ignore them. In the vast, vast majority of real-world cases (I'd rough guess at least 99.9%), the results from the higher value of Digits are closer to the truth. It may be a truth that you don't want to see, but it shows that the results from the lower value of Digits are completely unreliable. To ignore that is putting your head in the sand. That it makes the model fit the experimental data is never a sufficient reason. 

@mmcdara What the OP has written is correct when used in 2D Input. To verify this, copy the OP's input line V(r):= ... to a 1D Input (aka Maple Input) execution group. You'll get this (in Maple 2018):

V := proc (r) options operator, arrow, function_assign; piecewise(0 <= r and r <= R, -VV, R < r, 0) end proc

I think that the simulation that you want can be done in Maple itself. If you post details and formulas, we can get started.

@jefryyhalim I wouldn't say at this point that the behavior arises from the ODE. The red flag to me at the moment is this 1e-10, which I'll call epsilon.

To be clear about what we've learned so far: You know now that you can't use Digits = 10 with this small a value of epsilon, right? Digits:= 15 is a very good value to use because it's the largest value such that computations are done in the realm of hardware floats, which are numbers for which the very fastest computations are possible.

Now, you need to explain to me What is your purpose in using epsilon? I find it quite suspicious, and a likely cause of erratic, unstable, or undesirable behavior. I notice that you divide things by tan(Pi/2 +/- epsilon), which is essentially +/-infinity. This is exactly the kind of thing that is likely to cause Float(undefined) to appear.

@jefryyhalim Early in the worksheet, you set a variable to (essentially) evalf(Pi/2 + 1e-10). If Digits is 10 (its default value), then this is exactly equal to evalf(Pi/2). Surely that wasn't what you intended. I see that you had Digits:= 15, but you commented it out. Why? If you use Digits:= 15, do you still have the erratic behavior?