S:=a^4*b*c*x^2+a^4*b*c*y+a^4*b*c;

subs(t=a^4*b*c, collect(algsubs(a^4*b*c=t, S), t));

Use LinearAlgebra[DotProduct]  with the option  conjugate=false

Addition:  if you still use linalg[dotprod]  then use  'orthogonal'  option.

Example:

linalg[dotprod](vector([a,b,c]), vector([x,y,z]), 'orthogonal');

                                        a*x+b*y+c*z

Edited.

If a family of curves corresponds to the differential equation  F(x, y, y') = 0 , the orthogonal family will correspond to equation   F(x, y, -1/y') = 0 .

Example:

Equations  y'=2*x  and  -1/y'=2*x .  Their solutions  are  y=x^2+C  and  y=-1/2*ln(abs(x))+C .

Their plots:

A:=plot([seq(x^2+C, C=0..10)], x=-2..2, color=red):

B:=plot([seq(-1/2*ln(abs(x))+C, C=0..10)], x=-2..2, color=blue):

plots[display](A,B, scaling=constrained, view=[-2..2,-1..7]);

A more naive method:

e:= g^((2*(-sigma+k+1))/(-1+sigma))-tau^2:

ex:=(-sigma+k+1)/(-1+sigma):

subs(a=ex, factor(expand(eval(e, ex=a))));

Another way (without  ex):

e:= g^((2*(-sigma+k+1))/(-1+sigma))-tau^2:

factor(applyrule(g^(2*a::anything)=(g^a)^2, e));

Edited.

This is impossible, because in general, the roots of polynomials of degree higher than four can not be expressed explicitly in terms of their coefficients. Your problem can be solved numerically for any  w  by the procedure  Sol :

restart:

P:=w->w-10*z-1/5*z^2-1/200*z^3-1/500*z^4-z^5-1/1000*z^6:

Sol:=w->fsolve(P(w), z, complex):

Examples of use:

Sol(0);

Sol(1);

Sol(10);

If you only want the real roots, then remove  complex  option.

cat(2*x-5, ` > `, 3*x-6);

or even simplier:

 `2*x-5 > 3*x-6`;

For further changes, you can use  parse  command:

`2*x-5 > 3*x-6`;

parse(%);

solve(%);

Edited.

m := -t*(-1+t)^3/(t^3+2):

map(t->`if`(type(t,`^`) and op(2,t)>1, `$`(op(t)), t), [op(m)]);

                                   [-1, t, -1+t, -1+t, -1+t, 1/(t^3+2)]

Edited.

If the problem should be solved not a single time then the solution is useful to write down as a procedure. The parameter  P  is a set or a list of any number of polynomials. The procedure returns the list whose first element is the number of common terms of these polynomials, and the second one is the set of common terms.

CommonTerms:=proc(P::{set(polynom),list(polynom)})

local S;

S:=`intersect`(op(map(t->{op(t)},P)));

[nops(S), S];

end proc:

Example of use:

CommonTerms({x^2*y+x+y, x^3+x^2*y+y, y^2+x^2*y+y});

                                                [2, {y, x^2*y}]

In the code replace  Iteration  by  iteration  (or  iteration  by   Iteration). Maple distinguishes between lower and upper case.

Very strange bug! Use  Optimization[Maximize]  command as a workaround:

Optimization[Maximize](x^2+4, x = -1 .. 2);

                                            [8., [x = 2.]]

Carl showed how to find the general solution of a linear system. But if you want to modify a solution to the desired form, a special procedure is necessary. Here is my attempt :

LSSF:=proc(A::Matrix)

local Sol, Ind, n, N, Sol0, i, V;

uses LinearAlgebra;

Sol:=LinearSolve(A);

Ind:=indets(Sol);

n:=nops(Ind);

N:=op(A)[2]-1;

if n=0 then return Sol else

Sol0:=eval(Sol,Ind=~0);

for i from 1 to n do

V[i]:=eval(Sol-Sol0,{Ind[i]=1,op(Ind minus {Ind[i]})=~0});

od; fi;

print(x=Vector([seq(x[i],i=1..N)]),v[0]=Sol0,seq(v[k]=V[k], k=1..n));

x=v[0]+``(`+`(seq(C[k]*v[k],k=1..n)));

end proc:

Solution of the initial example:

R := Matrix([[1, -2, 2, 6, -6], [2, -3, 4, 9, -8]]);

LSSF(R);  # C[1] and C[2] are arbitrary parameters

 The procedure works for the arbitrary linear system which is given by an augmented matrix.

The absolutely exact (symbolic) value of  DthetaZero   can be gotten from the equation  Eq3  in my solution above:

DthetaZero:=eval((1/10)*sqrt(100*sin(theta(t))^2+1962*cos(theta(t))+1962), theta(t)=2*Pi/3);

evalf[20](%);

First we reduce the order of the equation. Using the conditions  theta(t)=Pi and diff(theta(t),t)=0  we find _C1 and then  DthetaZero . At this value there is instability when the pendulum reaches the top position. If the value  DthetaZero  is slightly less (DthetaZero1), the pendulum comes back, and if a little more (DthetaZero2), then the angle  theta  continues to increase (Plot1 and Plot 2):

restart;

Eq1:=sin(theta(t))*cos(theta(t))-9.8100*sin(theta(t))-(diff(theta(t), t, t)) = 0:

# The change diff(theta(t),t)=p. We get diff(theta(t), t, t)=p*diff(p(theta),theta)

Eq2:=p(theta)*diff(p(theta),theta)=1/2*sin(2*theta)-9.81*sin(theta);

sol0:=dsolve(Eq2)[2];

_C1:=solve(eval(sol0,[theta=Pi,p(theta)=0]));

Eq3:=eval(sol0,[p(theta)=diff(theta(t),t),theta=theta(t)]);

sol:=dsolve({Eq3,theta(0)=2*Pi/3},numeric):

DthetaZero1:=evalf[5](eval(rhs(Eq3),theta(t)=eval(theta(t),sol(0))));

DthetaZero2:=evalf[5](eval(rhs(Eq3),theta(t)=eval(theta(t),sol(0)))-10^(-7)+10^(-4));

Sol1:=dsolve({Eq1, theta(0) = 2*Pi*(1/3), D(theta)(0) = DthetaZero1}, numeric):

Sol2:=dsolve({Eq1, theta(0) = 2*Pi*(1/3), D(theta)(0) = DthetaZero2}, numeric):

plots[odeplot](Sol1,[[t,theta(t)-Pi],[t,diff(theta(t),t)]],t=0..14, color=[red,blue], axes=normal); # Plot 1

plots[odeplot](Sol2,[[t,theta(t)-Pi],[t,diff(theta(t),t)]],t=0..14, color=[red,blue], axes=normal); # Plot 2

For single use we can simply write

a:=0:  b:=1:  n:=4:

evalf[2]([seq(a+(b-a)/n*k, k=0..n)]);

                                                   [0., 0.25, 0.50, 0.75, 1.]

This also works for complex and symbolic cases.

restart;

[seq(a+(b-a)/4*k, k=0..4)];

                            [a, (3/4)*a+(1/4)*b, (1/2)*a+(1/2)*b, (1/4)*a+(3/4)*b, b]

a:=0: b:=1+I: n:=4:

evalf[2]([seq(a+(b-a)/n*k, k=0..n)]);

                                  [0., 0.25+0.25*I, 0.50+0.50*I, 0.75+0.75*I, 1.+1.*I]

If you can live with the extra parentheses, you can do so:

A:=<1/3,1/3; 1/3,1/3>;

1/3*``(A/~(1/3));