The coefficients of all the extraneous terms are relatively small. The relative error is ~ 1e-6 (just eyeballing it) over a small range of sigma2. However, it's true that an extra-degree term will ultimately cause unlimited error.  I'm working on some precise error computations, presented in a worksheet below.

The method used for the original determinant is unifloat. (So note that your guess that it was Gaussian elimination was incorrect.) The method used for Axel's method is fracfree. The details of each method that I summarize below can be viewed by setting infolevel[all]:= 5 before doing the computation.

Using the unifloat method, it first determines that the maximum possible degree of the determinant is 9. This estimate is very reasonable based on the maximum degree in every row and column. Both rowwise and columnwise estimates are 9 (as can be seen by inspection). Then it evaluates the matrix's floating-point determinant at 10 points and computes the interpolating polynomial. This method is nearly guaranteed to give you a dense degee-9 polynomial when done with floating-point arithmetic.

The fracfree method proceeds by relatively straightforward fraction-free Gaussian elimination. "Fraction-free" in this context means that no rational functions are used; it doesn't mean that integer coefficients are used.

Here's a worksheet showing the relative errors of the two methods.

Download poly_det_error.mw

plots:-implicitplot3d(
     (x-y)^2+(x-z)^2+(z-x)^2=3, x= -2..2, y= -2..2, z= -2..2,
     grid= [20$3], style= surface
);

Here is the solution. If you had Maple 2016, I would've used a DataFrame rather than a Matrix for the final output. That would make it eaier to read with column and row labels.

Download length_by_letter.mw

Another way, if you have trouble remembering the @:

evalf(Int(x-> h(g(x)), x= 1..2);

Using unapply with quotes on its first argument is pointless (and slightly wasteful). You could just as well define g via

g:= a-> fsolve(a*y^2 - sin(y), y=2);

So, I wonder where this post that recommends unapply is so that I can go correct it.

Here's another way:

plot([<x|y_1>, <x|y_2>], style= point, symbolsize= 24, symbol= [cross,box]);

To compute n-norm where n = 1 or 2, use

evalf(Int(abs(x(t - 2*Pi) - x(t))^n, t= a..b)^(1/n));

To compute the infinity norm, use

numapprox:-infnorm(x(t - 2*Pi) - x(t), t= a..b);

In either case, make sure to choose a and b appropriately. In the piecewise problem in the worksheet, for the infinity norm use

numapprox:-infnorm(eval(F, t= t- 2*Pi) - F, t= op([1,1,1], F)+2*Pi..op([-2,-1,-1], F));

     1.21380938300497*10^(-7)

Operands (op): Here's an abbreviated version of your piecewise:

F:=
piecewise(                                                                                                                         #op([0], F)
     5.442116637 <= t and t <= 8.035505284,                                                         #op([1], F)
     20.20246287+6.560150625*t-.5000000000*t^2-1.500000000*cos(t),   #op([2], F)
     8.035505284 <= t and t <= 8.583709290,                                                         #op([3], F)
     .4977821578                                                                                                               #op([4], F)

):
So op([1], F) is the conjunction 5.442116637 <= t and t <= 8.035505284. Negative numbers count back from the end so op([-2], F) = op([3], F) =  8.035505284 <= t and t <= 8.583709290. Now I explode the first conjunction to see what the suboperands are:

5.442116637 <= t     #op([1,1], F)
and                                #op([1,0], F)
t <= 8.035505284     #op([1,2], F)

Note that the main operator in any expression is its zeroeth operand regardless of whether it occurs at the beginning (as does piecewise) or in the middle (as does and). Now I explode the first inequality to see what its suboperands are:

5.442116637   #op([1,1,1], F)
<=                      #op([1,1,0], F)
t                          #op([1,1,2], F)

All the above information about operands applies to any (unevaluated) function, not just piecewise.

In the fourth equation, c is being used as a function applied to the expression in parentheses to its right. To avoid this, you need a space or an explicit multiplication operator (*) between the c and the following left parenthesis.

It's not an error, but your logp and logd are treated as independent constants. If you intend for these to be log(p) and log(d), then you need to enter them like that, with the parentheses.

Good Question. Vote up.

I recommend that you use command plot rather than odeplot for this purpose. The command plot uses an adaptive sampling scheme to avoid lacunarity, whereas odeplot is only controlling for numeric error (using dsolve's error-control parameters) rather than appearance. Here's an example:

randomize(7): #To select a specific example.
ODEsys:= {
     diff(y(t),t$2) = randpoly([sin(t), diff(y(t),[t$k]) $ k= 0..1], degree= 2),
     ((D@@k)(y)(0) = 'rand(-99..99)()') $ k= 0..1
}:
Sol:= dsolve(ODEsys, numeric, range= 0..2*Pi, maxfun= 10^5, output= listprocedure):
plot(eval([y(t), diff(y(t),t)], Sol), 0..2*Pi, numpoints= 10^3, adaptive= 6);
data:= [plottools:-getdata(%)]:
M||(1..2):= data[k][-1] $ k = 1..2;

Fear not: The solution procedures in Sol are still controlling for numeric error to the same extent that they are when using odeplot. The adaptive scheme used by plot means that at least numpoints points will be used. As a last resort, you can use plot's sample option to add some more points in trouble spots. It's very rare that you'd need to do that.

You have too many square brackets in your input of A. It should be like this:

A:=Matrix([[1/sqrt(2), 0],[0, 1/sqrt(2)]]);

or, even simpler,

A:= <1/sqrt(2), 0; 0, 1/sqrt(2)>;

The shoot package isn't a part of standard Maple. It's a third-party package that you need to find and download. Maple can usually solve BVPs by other means.

I told you not to transcribe the output or command prompts into your messages. Do you have a problem with that?

I'm not sure what you mean by "get", "answer", or by the 0 in your print command. You can "get" a plot by replacing the print command with

plot([seq(X2||k, k= 1..1)], 0..blt);

Is that what you want?

Note that X2 by itself means nothing; it always needs to be (X2 || something).

Unless you can correctly format the output, please don't copy the output to your posts. And skip the command prompts ("> ") also. It's easier for us to read and to use if you just give the raw input.

c = ln((1+sqrt(5))/2), exactly.

Proof: Let S[n] denote the set whose cardinality you've denoted as A[n]. Then S[n+1] = S[n] union {all members of S[n] that we're allowed to add n+1 to} = S[n] union (S[n-1] union~ {n+1}). The outer union is disjoint because no member of S[n] contains n+1. Thus A[n+1] = A[n] + A[n-1]. It is well known that the asymptotic ratio of successive terms of this recursive sequence is (1+sqrt(5))/2 regardless of initial conditions, but I can give a more-detailed proof if you want.

P.S.: Some Maple code:

S[-1]:= {{}}:  S[0]:= {{}}:
for n to 22 do
     S[n]:= S[n-1] union map(`union`, S[n-2], {n})
od:
seq(nops(S[n]), n= 0..22);

     1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368

Maple won't allow the union~ that I used in the exposition because ~ can't be used between container objects of different sizes; so I needed to use map(`union`, ...) instead.

I can't verify that this will work without seeing your function. Let's say that the function is f(x). Try

simplify(expand(f(x+2*Pi) - f(x)));

If you get back 0, you can consider that to be a verification that the function is 2*Pi periodic. If you don't get back 0, try making some assumptions on x, such as 

simplify(expand(f(x+2*Pi) - f(x))) assuming x >= a, x <= c;

where c is some critical point between a and b, then

simplify(expand(f(x+2*Pi) - f(x))) assuming x >= c, x <= b;

I don't know how to address your more-general question: How to find a period of a function.

First, don't use . for multiplication; use * instead.

N:= 4;
y:= sum(A[2*n]*cos(2*n*x), n= 0..N);
eq1:= diff(y,x,x)+(a+2*q*cos(2*x))*y;
eq2:= map(combine, eq1, trig);
eq4:= map(collect, [coeffs(eq2, [seq(cos(2*n*x), n= 1..N+1)])], [seq(A[2*n], n= 0..N)]) =~ 0;

     [a*A[0]+q*A[2] = 0, (a-4)*A[2]+2*q*A[0]+q*A[4] = 0, q*A[2]+(a-16)*A[4]+q*A[6] = 0,
     q*A[4]+(a-36)*A[6]+q*A[8] = 0, q*A[6]+(a-64)*A[8] = 0, q*A[8] = 0]

Note that there's a cos(10*x) term and that the coefficient of 1 is always returned by coeffs.

I haven't tested the crashing aspect yet because obviously your int should be Int, and (not so obviously) the eval(...) right before the end do should be evalf(eval(...)). Make those changes and you'll get numeric answers immediately.

I don't see what thread safety has to do with it. This is not parallel processing code.