Your worksheet is quite complicated, so maybe you can give some more information about the matrix or provide a simpler example. In particular, are all the entries numeric? is it sparse? Without that information, I can make only the general suggestion

You want FF3 = FF2^(-1).FF1

In general, finding the matrix inverse is bad practice here - you should just solve FF2.FF3=FF1 for the unknown matrix FF3.

Just call LinearSolve with A=FF2 and B=FF1:

FF3 := LinearAlgebra:-LinearSolve(FF2,FF1);

Things will be slow for Digits=30. If you can use Digits less than evalhf(Digits) (usually 15), then the numerical computations will be done in hardware and will be much faster.

Since a=a means any a is acceptable, and you want a>0, then a=a,b=0,c=0 and d=d is not a solution. But any d is OK, so I think a>0,b=0,c=0,d=d is a missing solution.

Download ineq.mw

There is only one value for each element. According to ?ScientificConstants,properties, "for those elements with several solid forms, the density is of the most stable crystal form". This is the usual convention for thermodynamic tables, where the most stable allotrope is the standard form (one with the enthalpy of formation defined as zero). The exception is phosphorus, where white phosphorus is the standard form, even though red phosphorus is more stable (because the red phosphorus structure is not well defined). In Maple too, the phosphorus density given is that for white phosphorus (in contrast to the help statement).

nops of a Vector is always 3 - you want numelems here.

As in the help page, Maple outputs the normal form without the D. As per the Wikipedia page, the D can be found from L by taking the squares of the diagonal entries. Following the prescription there, if L is the matrix from Maple's output:

S := DiagonalMatrix([seq(L[i, i], i = 1 .. Dimension(L)[1])]);

D:=S^2;

New L:

Lnew:=L.S^(-1);

Following Carl's suggestion, the algorithm seems to give eigenvectors with the symmetry you want, except for those corresponding to the pair of zero eigenvalues, and the matrix of eigenvectors produced is orthogonal (as expected for a symmetric matrix). Furthermore the eigenvectors corresponding to the zero eigenvalues both seem to have their 2nd half entries the negative of the first half entries (see below for an example), say [x,-x] and [y,-y] (this is sloppy notation for column vectors stacked above each other).

Now the only possibility is that there are linear combinations of the eigenvectors for the zero eigenvalues that have the symmetry you want. Now any linear combination of [x,-x] and [y,-y] is [z,-z] and so any two orthogonal linear combinations must be [z,-z] and [z2,-z2]. But you want these [v1,v2] and [v2,v1] which can only occur if v2=-v1 so [v1,-v1] and [-v1,v1]. But these are multiples of each other and so not linearly independent.

Perhaps I'm missing something, but are you sure that this is possible?

Download SymmetryEigenvectors2.mw

Edit: Note that NullSpace(M) also gives eigenvectors of the form [x,-x], so the above isn't an artifact related to the fact that the two eigenvalues are small but not exactly the same.

Use Pi not pi for the mathematical constant 3.1416...; pi is just a symbol without a value. The Greek palette gives Pi not pi. 

Move from Text to Math and then choose 2D Math (not 2D Math input), as shown above.

@Rouben Rostamian's answer has the details. Explore allows you to easily vary y and see what happens to the plot. If you do this to the derivative, you see where the zeroes are:

Explore(plot(diff(x^2+y*(1-x)^(3/4),x),x=0..1),y=0..2.0);

or you can look at the plot and its derivative together:

f:=x^2+y*(1-x)^(3/4):Explore(plots:-display(<plot(f,x=0..1)|plot(diff(f,x),x=0..1)>),y=0..2.0);
 

firint in the DEtools package can calculate first integrals, which may be helpful. For your second order ode, it successfully produces a first order ode.

As I indicated to the OP on the other thread, if fsolve doesn't find an answer, it returns NULL (as in the response by @tomleslie) or unevaluated (as in the response by @mmcdara). If you want to test for both these possibilities, then you want

ans:=fsolve(x=x+1,x);

if ans=NULL or op(0,ans)=`fsolve` then "no" end if;

Note this handles also the multivariable case, where a valid solution won't just be a number or sequence of numbers, but something like {x=1, y=3}.

Testing NULL output from fsolve should work here, though sometimes fsolve just returns unevaluated.

Seems to me it is working, just fsolve is really slow. You can speed it up by using the last solution as the initial estimate (but convergence of your outer loop is still really slow).

Try

Q[h+1] := fsolve(eval(f[3], {R = R[h], S = S[h]}) = 0, Q = Q[h]);
if Q[h+1] = NULL then print("breakQ"); break end if;

If you work out Theta[k] and Phi[k] for k=2 nd 3, then you can use a single loop from 4 to M for Theta[k], Phi[k], F[k]. (I shouldn't really cut and paste to duplicate lines here.)

Solution_1-dharr.mw

as @mmcdara pointed out, you van specify ranges in which to search for solutions. You can also use the avoid option to not get the same solution, e.g. after ans1:=fsolve({f, g}); you can use fsolve({f, g},{x,y},avoid=ans1); to get the second solution.

plots:-complexplot(exp(I*theta),theta=0..2*Pi);