Your procedure can be further improved like this:

sim:= proc(                                                             (*Notes*)
   x::posint,                                                           (*1*)
   n::And(posint, satisfies(n-> n >= x)),                               (*2*) 
   o::identical(exactly, minimum, maximum):= 'exactly'
)
uses RT= RandomTools;                                                   (*3*)
local
   counter, roll,
   relation:= Record('exactly'= `=`, 'minimum'= `>=`, 'maximum'= `<=`)  (*4*)
;
   for counter do                                                       (*5*)
      roll:= RT:-Generate('list'('integer'('range'= 1..6), x))          (*6*)
   until relation[o](add(roll), n);                                     (*7*)
   roll, counter                                                        (*8*)
end proc:

By far the most important of these improvements is (*8*) that I removed print. A procedure should never use print to return its value. You may use print to return supplementary information. If a procedure tries to return its value with print, then you won't be able to assign that return value to a variable, and you won't be able to use the procedure effectively inside other code.

Other notes:

(*1*) It's usually a good idea to put a type restriction on all parameters (posint means positive integer).

(*2*) That type restriction can include a dependency on other parameters. Here, n is restricted to be a positive integer greater than or equal to x. This can usually be done without using the depends modifier, whose documentation is quite weak (so I recommend that you completely avoid that while you learn the rest). The type satisfies is used to turn any boolean condition into a type.

(*3*) It's my personal preference that package names not be declared 

uses package name;

but rather

uses abbreviation= package name;

This makes it absolutely clear which package, if any, the members come from.

(*4*) Locals can be assigned values directly in their declaration.

A Record is one of Maple's many indexable container data structures. It can be used when the indices are all symbols, as are exactly, minimum, and maximum. The target values can be anything.

When operators (such as <>, <, >) are referenced without operands (or when used in prefix form), they need to be put in back quotes: `<>`, `<`, `<`. Don't confuse these with the forward quotes used to prevent evaluation, as in 'exactly'.

(*5*) When a counter needs to be incremented for every pass through a loop, you should use for to do it. You can include a while clause on the for line, but in this particular case, I thought that an until clause was better.

(*6*) Note that RandomTools:-Generate can generate an entire list in a single call by using this syntax. The second argument to list(...) is the length of the list, x in this case. This is more efficient than using repeated calls to Generate.

(*7*) The until is a clause of the do statement (as in for ... do). It's not a statement that can stand alone. It replaces the end do or od. It's like the opposite of while, but unlike while, it forces the loop to be executed at least once.

The relational operator is being used here in prefix form.

When add (or seq or mul) are used on a list or Vector, there's no need to use an index.

@digerdiga Your A and B return lists of different lengths when given the same input?? I don't understand the point then. You'd have the same problem if they were specifed as linear.

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