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);

put the plot in a variable, e.g., p1:=plot(x^2,x=0..3):

and then lprint(p1);

Even though the first two seem to be interpreted as 2-D integrals, the correct syntax for 2D integrals needs the variables and ranges in a list (or nested 1D integrals), like

simplify(Int(integrand, [phi = 0 .. 2*Pi, r = 0 .. 2]));

and then you get the correct answers. 

You have two first order equations, so if you know omega[0] and sigma, you need only two initial conditions. Put them as a sequence not as a list to solve the syntax error. So give omega[0] and sigma values then

ICS:=C(0)=0,A(0)=1:
sol:=dsolve({C_t,A_t,ICS},{C(t),A(t)},type=numeric):
pRange:=0..20:
plots:-odeplot( sol, [t, A(t)],t=pRange, numpoints=10000 );  #note A(t) not A

If you want to find the values of omega[0] and sigma by specifying two additional conditions, then I think that only works with boundary value problems, and you have only initial conditions. Note also the syntax for derivative conditions

is D(C)(0)=0 for derivative zero at time zero 

I think dchange existed in Maple V, so try

integral := Int(2*(sin(theta)/cos(theta))^(2*p-1), theta = 0 .. (1/2)*Pi);
PDEtools[dchange]({theta=arctan(t)},integral);
value(%);

gives Pi*csc(Pi*p) (in Maple 2017)