Not only is it not a function, but it is not a PDF (the F in PDF stands for function).  This is what's called a "distribution of mixed type"  (partly discrete, partly continuous).   See e.g. www.math.uah.edu/stat/dist/Mixed.xhtml

When you give your equation to solve (without specifying the value of Z) you get a single expression involving RootOf a polynomial.  Presumably many of the roots are complex, so even if you use allvalues to get the individual roots max won't do well with this.  Somehow (I'm not sure why), if you give Z a numerical value solve takes a different route and returns numerical values.

When you give your equation to solve (without specifying the value of Z) you get a single expression involving RootOf a polynomial.  Presumably many of the roots are complex, so even if you use allvalues to get the individual roots max won't do well with this.  Somehow (I'm not sure why), if you give Z a numerical value solve takes a different route and returns numerical values.

Yes, the subspace associated to any isolated cluster of eigenvalues should be OK: the eigenprojection on it could be obtained by a contour integral
1/(2*Pi*I)*int_C (z*IdentityMatrix - A)^(-1) dz 
where C is a contour enclosing those eigenvalues and no others, and this is an analytic function of the entries of A as long as the eigenvalues stay away from C. 
In principle Maple could compute the eigenprojection using LinearAlgebra:-MatrixFunction, but I don't know how numerically stable that would be.

Yes, the subspace associated to any isolated cluster of eigenvalues should be OK: the eigenprojection on it could be obtained by a contour integral
1/(2*Pi*I)*int_C (z*IdentityMatrix - A)^(-1) dz 
where C is a contour enclosing those eigenvalues and no others, and this is an analytic function of the entries of A as long as the eigenvalues stay away from C. 
In principle Maple could compute the eigenprojection using LinearAlgebra:-MatrixFunction, but I don't know how numerically stable that would be.

The case of close but not equal eigenvalues can be the source of very severe numerical instability.  Consider the following two matrices, which are very close together when epsilon is very small.

A:= << 1 -epsilon , epsilon > | <epsilon,1 - epsilon>>;
B:= << 1 - 2*epsilon, 0> | <0, 1>>;

Both have eigenvalues 1 and 1-2*epsilon which are very close together, but the corresponding eigenvectors of A are
<1,1> and <-1,1> while those of B are <0,1> and <1,0>.

The case of close but not equal eigenvalues can be the source of very severe numerical instability.  Consider the following two matrices, which are very close together when epsilon is very small.

A:= << 1 -epsilon , epsilon > | <epsilon,1 - epsilon>>;
B:= << 1 - 2*epsilon, 0> | <0, 1>>;

Both have eigenvalues 1 and 1-2*epsilon which are very close together, but the corresponding eigenvectors of A are
<1,1> and <-1,1> while those of B are <0,1> and <1,0>.

You're solving P(s) = 2.5, not P(s) = .25.  The answer Maple returned was quite correct: that equation has two complex solutions, but no real solutions.

You're solving P(s) = 2.5, not P(s) = .25.  The answer Maple returned was quite correct: that equation has two complex solutions, but no real solutions.

Simpler, I think, is to use eval.

> eval ( [temp, x1, x2, y1, y2], sol);

[234.9000000, .7633930119e-14, .2397737599, 1402.995860, -47636.68497]

or if you want the list of equations with variable names

> [seq(t = eval(t, sol), t = [temp, x1, x2, y1, y2])];

[temp = 234.9000000, x1 = .7633930119e-14, x2 = .2397737599, y1 = 1402.995860, y2 = -47636.68497]

Simpler, I think, is to use eval.

> eval ( [temp, x1, x2, y1, y2], sol);

[234.9000000, .7633930119e-14, .2397737599, 1402.995860, -47636.68497]

or if you want the list of equations with variable names

> [seq(t = eval(t, sol), t = [temp, x1, x2, y1, y2])];

[temp = 234.9000000, x1 = .7633930119e-14, x2 = .2397737599, y1 = 1402.995860, y2 = -47636.68497]

> subs(A=Matrix([[1,2,3],[4,5,6],[7,8,9]]), X = <x1,x2,x3>, B = <b1,b2,b3>, A.X = B);

> subs(A=Matrix([[1,2,3],[4,5,6],[7,8,9]]), X = <x1,x2,x3>, B = <b1,b2,b3>, A.X = B);

However, in many cases it does not work.  For example, it won't work with a Matrix.

RTABLE(73297976,MATRIX([[1, 2], [3, 4]]),Matrix)
 

And here's another example that gets messed up:

1/x^2+2/x+3

It's not at all clear what you mean.  You have one equation with three variables: A, B and x (since you've said H=800). 

Now you say you want x = 0 .. pi/2, but  there are no constant values of A and B that will make this equation true for all x in that interval (or any other interval of nonzero length), because sin(x), sin(2*x) and sin(4*x) are linearly independent functions.  You could solve for A as a function of B and x, or B as a function of A and x.