Although pasting works sometimes, your learning is much improved by explicitly typing examples even when pasting works. Perhaps as you typed, you'd realize that what's pasted above is nonsense because it's missing both the multiplication and exponentiation operators.

Obtaining non-real floats (i.e., with nonzero imaginary parts) when you expect only reals is not a sign of an error or a cause for concern. There are two reasons that you might get these: round-off error or ignoring symbolic constants.

Round-off error: The decimal approximation of some real constants often necessarily involves detours into the complex numbers.  For example, ((-1)^(1/3))^3 is clearly -1, but if I include a decimal point:

((-1)^(1/3.))^3;
      -.9999999999+3.865045973*10^(-10)*I

The small magnitude of that imaginary part is a sign that that imaginary part is simply due to round-off error and can be ignored. My rule of thumb for this is if you're expecting a real number and

Im(z) <> 0 and Re(z) <> 0 and ilog10(Im(z)) < 3-Digits and ilog10(Im(z)) - ilog10(Re(z)) < 3-Digits

and these relations remain true when you increase Digits, then the imaginary part can be ignored and/or discarded. (It's just a rule of thumb: there are a few situations where it doesn't work.) There are easy ways to remove these imaginary parts which we can go over when the need arises. At this developmental phase of your coding project, you should acknowledge their presence and mentally ignore them without explicitly removing them.

Ignoring symbolic constants: I don't know where you're at in your mathematical education, but you're probably at least vaguely aware of the following (each of these is either covered in a precalculus course  or a first course in differential equations).

1. The general solution of a homogenoous linear ODE with constant coefficients involves finding the roots of a univariate polynomial with the same coefficients and with the degree of each term being the differential order of the corresponding term in the ODE.
2. That if the roots r[k] of that polynomial are distinct then the general solution is y(t) = Sum(C[k]*exp(r[k]*t), k= 1..n), where n is the differential order of the ODE and the C[k] are constants which can't be specified until initial or boundary conditions are provided. (Only minor adjustments to that are needed when the roots aren't distinct.)
3. The roots of a polynomial are in general non-real even if the coefficients are real. 
4. If the coefficients are real, then any non-real roots occur in complex-conjugate pairs (e.g., x+I*y and x-I*y).
5. For real t, exp(I*t) = cos(t) + I*sin(t). (That's still true if t is non-real, but we only need to consider the real case here.)
6. The sum of any two solutions of a linear homogenous ODE is also a solution.
7. There is no formula or algorithm for finding the exact roots of a general polynomial of degree greater than 4. (This is one reason why you often see polynomial roots expressed in Maple's RootOf form.) Indeed, it's been proven that no such formula or algorithm will ever be found. There are algorithms to find decimal approximations of the roots to arbitrary accuracy.

If you put these facts together, the result is that the real solution to a real-coefficient ODE without initial/boundary values might be expressed as y(t) = sum(C[k]*exp(r[k]*t), k= 1..n) where some of the r[k] are manifestly non-real (not of the "round-off error" variety described above), and in complex-conjugate pairs. When the values of the C[k] are determined, some of them will also be non-real, but they will be such that the overall solution is real valued, upto the limitations of round-off error.

Numeric BVP solver: I think that your problem can be solved by Maple's numeric BVP solver, which is accessed by 

Sol:= dsolve(
   {some ODEs,
    some boundary conditions (BCs) expressed at the two endpoints of a closed real interval of the independent variable
   },
   numeric,
   possibly some other numeric options (as needed)
);

You should at least try. Symbolic parameters can be included as long as you have exactly 1 BC for every symbolic parameter plus the total differential order of the ODEs. The BCs can be quite complicated as long as the only explicit values of the independent variable that are used are the two endpoints of the interval. You can use a parameter for the eigenvalue. The solution of the system will include real numeric values for all parameters. This technique completely avoids non-real values. See help page ?dsolve,numeric,bvp.

Lengthy output messages: The message "Length of output exceeds ..." is not an error message nor even an indication that there was any snag in computing the solution that you requested. The message is merely there to warn you that the computed output is so long that you may not want it displayed on your screen, or at least you may not want it in prettyprinted format. The easiest way to look at the solution is

lprint(%);

which'll give it in plaintext form. However, this is always a good time to ask yourself "Is it reallly worth looking at? Can I analyze it in some other way (numeric evaluation, plot, breaking it into parts, etc.)?"

For x close to 0, x is a close approximation of sin(x). You can see that from its Maclaurin series. So all you're seeing are shifted approximations of sqrt(10).

My workaround is similar to Tom's, and additionally it'll continue to work if and when Maplesoft changes the parameters, regardless of how they're changed (even if it's more complicated than just switching their order):

MyGammaD:= proc(k, theta)
uses St= Statistics;
local a,b, GD:= 'GammaDistribution'(a,b);
   eval(GD, solve({St:-Mean(GD) = k*theta, St:-Variance(GD) = k*theta^2}, {a,b}))
end proc:

I'm not claiming to know whether Maplesoft has ever changed the parameters. The above works regardless of changes.

Even more bothersome to me than this thing with GammaDistribution is that there doesn't seem to be standardization across mathematical software and literature as to whether the second parameter of Normal is sigma or sigma^2.

I don't know what you're trying to achieve with the unevaluation quotes, but they have no effect whatsoever when used on a procedure parameter. Other than that, it looks like you're trying to match expressions that are invariant under op. Any conditional can be made a type by using type satisfies. So, the type is satifies(x-> x = op(x)). This won't error out if op(x) is multi-term, so there's no need for square brackets.

With tools like unapply and ->, there's no practical limit to the level of "meta differential operators" that you can create, each capable of being iterated with @@. In this example, I define an operator Newton which takes an "operator template" (f,c) and returns the differential operator that constructs the Newton's method iteration operator for solving f(x) = c. This final operator can be iterated with @@. Each level of these operators is only one line of code.

Newton:= proc(f,c) local x; subs(_c= c, _c-> unapply(x - f(x,_c)/D[1](f)(x), x)) end proc:
Sqrt:= Newton(((x,c)-> x^2-c), c): 
(Sqrt(2)@@4)(1.5); #1.5 is arbitrary starting point.
                          1.414213562
%^2;
                          1.999999999

The output produced by Newton itself is very ugly to look at. But the result of Sqrt(2) is what you'd expect.

Regardless of whether its number of elements is even or odd, the median^[1] of a nonempty^[2] sorted list X equals

add(X[[floor,ceil]((nops(X)+1)/2)])/2

That idea is used in the following, where I also address an issue that you'll encounter when your input has symbolic real constants that you may not want to clobber with evalf in the output:

#Note that the default sort order doesn't respect the usual order
#of real numbers when the numbers are symbolic:
sin~([$1..4]): (sort, evalf)(%);
     [sin(1), sin(2), sin(3), sin(4)], [0.8414709848, 0.9092974268, 0.1411200081, -0.7568024953]

#We can get the desired order by sorting with evalf:
(x::list)-> ((S,n)-> (S, add(S[[floor,ceil]((n+1)/2)])/2, add(S)/n))(sort(evalf(x)), nops(x)):
   
%(sin~([$1..4]));
     [-0.7568024953, 0.1411200081, 0.8414709848, 0.9092974268], 0.4912954964, 0.2837714810

#We can do the same thing plus not clobber the input
#by using sort(..., output= permutation):
(x::list)-> 
   ((P,n)-> (x[P], add(x[P[[floor,ceil]((n+1)/2)]])/2, add(x)/n))
      (sort(evalf(x), 'output'= 'permutation'), nops(x))
:   
%(sin~([1,2,3,4]));
                                      1          1         
    [sin(4), sin(3), sin(1), sin(2)], - sin(3) + - sin(1), 
                                      2          2         

      1          1          1          1       
      - sin(1) + - sin(2) + - sin(3) + - sin(4)
      4          4          4          4 
      
evalf([%]);
     [[-0.7568024953, 0.1411200081, 0.8414709848, 0.9092974268], 0.4912954964, 0.2837714811]

^[1]I mean median in the sense that it's taught in elementary math, not the more-technical sense discussed by @sand15 .

^[2]And do we really want to give serious consideration to defining the median of an empty list?

The vertical axis's tickmarks are positioned at -1/2, -3/2, ..., -(2*n-1)/2. Thus, this works:

restart:
n:= 10:
RM:= LinearAlgebra:-RandomMatrix(n):
NewTickMarks:= [($1..n)+~1/2 =~ A||(1..n)]:
Statistics:-HeatMap(RM, tickmarks= [NewTickMarks, (1-lhs=rhs)~(NewTickMarks)]);

You have 4 real-valued dimensions: x1, x2, x3, and f1. You're only allowed 3. If you set one of those 4 to a constant value, then you could plot the remaining 3 as a level surface. Then you could use time as a 4th dimension and make an animation. Sometimes I use gradations of color to represent extra dimensions.

On the left side of the assignment in the double loop, you have uppercase P. I think that you meant lowercase p.

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.