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

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