Unless for polynomial systems for which you can specify the "maxsols" number of solutions to find, fsolve returns only one.
And of course the solution you get depends on your starting point (see fsolve[details] help page) or on the search interval you set (same help page).

An alternative is to use the function "Isolate" from package RootFinding.

You will find a few results in the attached file.

RootFinding.mw

Not knowing all the subtitlies of the way Maple works, my first reflex when the arrow operator doesn't return the expected result is to use "unapply".
It goes well here:

for i from 1 to 3 do
   f[i]:=unapply(P(i,1,1,x), x);
end do;


To Carl, vv and others: what are the pros and cons of all the solutions proposed here ?

Answer "by hand"

product(1-1/k, k = 2 .. infinity)     
= limit(product(1-1/k, k = 2 .. N), N=+infinity)    
= limit((1-1/2)*(1-1/3)*...*(1-1/N)), N=+infinity) 
= limit((1/2)*(2/3)*(3/4)*...*(1-1/N)), N=+infinity)   # the denominator of the term n is the numerator of term n+1
= limit(1/N, N=+infinity) 
=0

Simpler than expanding (1-1/2)*(1-1/3)*.... = 1 - 1/2 - 1/3 - 1/4 - ... + 1/2*1/3 + 1/2*1/4 + ... + 1/3*1/4 + ... + triple products + quadruple products + and so forth...



                 

Before receiving a more astute answer here is a solution

Matrix.mw

I'm not sure to understand what you mean by "there is a spur at node 5"?

Mor of this, with your  [[14, 16, 1], [1, 2, 2], [2, 3, 7], [3, 4, 1], [4, 5, 3], [5, 6, 6], [6, 7, 4], [7, 1, 3], [8, 9, 7], [9, 10, 3], [10, 11, 3], [11, 12, 1], [12, 13, 6], [13, 14, 3], [14, 8, 3], [5, 10, 3], [7, 8, 3]] definition there are 15 vertices and 17 edges, not 14 and 16 (if I'm not mistaken, unless you count vertex 16 as a particular "source node" ?).

Besides, when you write [1, 2, 2], for instance, does it mean "there exist an arc from 1 to 2 of cost 2" or "there exist on arc from 1 to 2 and one arc from 2 to 1, both of cost 2"?

I guess you already know that GraphTheory contains a DijkstrasAlgorithm procedure (maybe it doen't suit you because of this 
" spur at node 5" you mentioned). 
Just in case I send you this little worksheet

Dijsktras.mw 

Maybe writedata could be useful (it depends on the formats uou use)

Here is an example (to write in a file just change "terminal" by the name if the file)
N := 5:
A := Vector[row](N, i->i);
writedata(terminal, [convert(A, list)], [seq(float, n=1..N)]);
1    2    3    4    5


fd := fopen(MyFile):
writedata(fd, [convert(A, list)], [seq(float, n=1..N)]):
close(fd):

Hint: The product will be complicated algebraically and you will need to apply first expand, then simplify to reduce the product of the eigenvalues to an integer. First try to do for a single matrix , then make a loop to do it 8 times.] 

If this is the hint your teacher has given you, maybe it's not necessary to mention it here: let the people think they do not do you homework.
However this is a feeble hint: if you use RandomMatrix the resulting matrix contains only numbers (integers or floats depending on the  way you describe the generator [see the RandomMatriwx help page]).
Then Eigenvalues will return 3 numbers and their product trivial (no need to expand anything).
restart:
with(LinearAlgebra):
A := RandomMatrix(3, 3, generator=0..1.0):
EigA := Eigenvalues(A);
EigProduct := mul(Re~(convert(EigA, list)));
Determinant(A)

PS: if you want to verify that the product of the eigenvalues equals the determinant, it would be beter to generate a symbolic matrix :
A := Matrix(3, 3, symbol=a):

Question 1:Are you looking for a solution for this specific polynomial or for a solution metho for more general polynomials ?
Question 2: Whjat is for you the shortest monomial coefficient among these two: p+q+r and p*q*r ?
Question 3: what about if two monomoal coefficients have the same "length" according to your definition ?

For your specific case tou will find an answer in the attached file

ShortestMonomialCoefficient.mw

Just replace 
( Vector[row]
([seq
by
(Matrix(
(3, 2, [seq

I have no solution to propose, only an experience to share.

From my experience cmaple (Windows) and maple (unix/mac OS) have some limitations.
The use of ExcelTools seems effectively to be one of them. 

Some features contained in "recent" packages are not handled too :

• Using functions from Grid package generate errors (Maple 2015, Maple 2016 with Windows)
• Some functions from DocumentTools are not supported: this is rather normal when they refer to displays and plots, but RunWorksheet doesn't work.
• If I'm not mistaken plotdevice  doesn't work too (at least in the situations I used it).

It's rather difficult to find a list of functions (if it even exists) that do not give the expected result when used with cmaple/maple, so, unfortunately, you eventually discover that by yourself when cmaple fails.

I did my grieving of all this: my opinion is that cmaple/maple are historical commands the evolution of which doesn't follow the evolution of xmaple. So I have adopted this pragmatic strategy: replace all the functions that are "not supported" by others which are.

For instance, when I faced a problem with ExcelTools:-Export, I replaced this command by a writedata/write.ImporttMatrix one (depending on the data to export).
Maybe you can do the same for ExcelTools:-Import: open you file with Excel, save it as a text file, and read it with a read/readdata ... command.

I hope someone here will give you better advices than I did
 

You will find in the attached file something that might help you or give you hints to go further.
I do not pretend it's the most elegant way to proceed.
I did not use your worksheet but a notional example gifen in the CodeTools:-Usage help page

CodeToolsUsage.mw

dsolve returns the exact solution but it appears to be definitively untractable.

I propose you the solution obtained after replacing the integer coefficients by floats


MaybeThis.mw

restart:
sqrt(sqrt(p^2+1)-1)*sqrt(sqrt(p^2+1)+1):
combine(%, radical)  assuming p > 0:
simplify(%)   assuming p > 0

restart:
sqrt(sqrt(p^2+1)-1)*sqrt(sqrt(p^2+1)+1):
combine(%, radical)  assuming p < 0:
simplify(%)   assuming p < 0

Maybe this will suit you if the term "tensor" stands for "a 4 dimensional table"

H := Array(seq(1..3, k=1..4)):
H[1, 1, 1, 1] := 10:
H[1, 2, 2, 1] := 20:

But if you want to use "true" tensors you can have a look to the Physics or DifferentialGeometry packages

Please see my other replies.
I've posted 
dsolve({ic, odes, odex}, type = numeric, range = 0 .. 1, events = [[s(tau)-1, halt]])

instead of 
dsolve({ic, odes, odex}, type = numeric, range = 0 .. 1, events = [[s(tau)-const, halt]])

APOLOGIES



About your original question
 [ [ [s(t), a*arcsinh(2/a) < s(t)], halt] ] means : halt dsolve when s(t) = 0 AND a*arcsinh(2/a) < s(t)]

I don't know what is the error Maple returned, but there could be potential conflict between these two conditions.

Is "a" a numerical value ?
If not it can't work until you have used the option parameters (see my first reply)
... or given "a" some value...
Look to attached file after thet "EXTRA MATERIAL" title

a := 3.0:
dsolve({ic, odes, odex}, type = numeric, range = 0 .. 1, events = [ [ [s(tau), s(tau) > a*arcsinh(2/a)] , halt]] )

dsolve.mw