A more detailed investigation shows that  f  has not the values  -1 or 1 anywhere.  Therefore, the range  of this function is the open interval  -1 .. 1 :

f := sqrt(x^2+x+1)-sqrt(x^2-x+1):

solve(diff(f,x)>0);  # The function is strictly increasing on the whole axis  x

limit(f, x=-infinity);

limit(f, x=infinity);

minimize(f, location); 

maximize(f, location);

plot(f, x=-10..10);

Joe, your first method is very witty, vote up! Here is another simple method using ListTools  package tools. It does exactly what the OP wants to:

L := op~([g[ ]]);
S:={L[ ]};
for x in S do print(cat(x,`  appear  `, nops([ListTools[SearchAll](x, L)]))); end do;
 

A:=2*Pi^2*(-1+delta)/r^(1+delta);
-evalindets(A, 'And(`+`, satisfies(f->sign(op(1,f))=-1))', t->-t);

Tom, your procedure will not work if one of the ends of the range is equal to infinity. Here's my version:

ConvertTo_range:=proc(S::{set,list})

local T, T1;

T:=map(t->[op(t)],S);

if nops(T)=2 then

T1:=seq(op(remove(s->type(s,name),T[i])), i=1..2);

return min(T1)..max(T1) else

if type(op(T)[1],name) then return -infinity..op(T)[2] else op(T)[1]..infinity fi; fi;

end proc:

Example of use:

L:=[{ arccos(3/4) < theta},

{arccos(-3/4+(1/4)*sqrt(13+16*sqrt(2)))+Pi <= theta, theta < 2*Pi-arccos(3/4)},

{theta <= 3*Pi-arccos(-3/4+(1/4)*sqrt(13+16*sqrt(2))), 2*Pi+arccos(3/4) < theta}]:

ConvertTo_range~(L);

evalf(%);

Example:

select(s->s[2]="e" and s[4]="d", ["trend", "fends", "noods", "sends"]);

x:=<x1(t), x2(t)>;

A:=<1,1; 1,-1>;

sys:=diff~(x,t)=~A.x;

ic:=eval(x,t=0)=~<1, -1>;

dsolve(op~(convert~({sys,ic}, list)));

or using DEtools (shorter)

x:=<x1(t), x2(t)>;

A:=<1,1; 1,-1>;

DEtools[matrixDE](A, t, method = matrixexp);
x=%[1].<1,-1>;

Another workaround in 2d input  to use  seq  command:

A := Array(seq(1 .. 3, i = 1 .. 3));

P  procedure solves the problem for 3 equations with 3 variables  x, y, z:

P:=proc(f, g, h)
if (not depends(f, z) and not depends(g, z)) or (not depends(f, y) and not depends(h, y)) or (not depends(g, x) and not depends(h, x)) then false else true fi;
end proc:

Example of use:

P(x*y*theta[5]+x*theta[2], x*y*theta[15], x*theta[22]+y*theta[23]+z*theta[24]);

                                                            false

Addition. If you want the solution for an arbitrary number of variables and equations, it is necessary  more precisely determine that the procedure must return in the different cases.

To use  NumberTheory[RepeatingDecimal]  command a decimal must be specified as a fraction. The procedure  ToFraction  converts any decimal with a recurring part to a rational number. You can then use  RepeatingDecimal  command.

Procedure parameters: a - part of a recurring decimal before of the repeating part, R  is a repeating part defined as a string.

ToFraction := proc(a::float, R::string)

local m, n, b;

m := SFloatExponent(a);

n := length(R);

b := add(parse(R[i])*10^(n-i), i = 1 .. n);

convert(a, fraction)+sign(a)*10^(m-n)*(sum(b*10^(-n*(k-1)), k = 1 .. infinity));

end proc:

Examples of use:

ToFraction(0., "142857"); 
NumberTheory[RepeatingDecimal](%);

ToFraction(1.6, "12"); 
NumberTheory[RepeatingDecimal](%);

ToFraction(1., "0102"); 
NumberTheory[RepeatingDecimal](%);

Put braces or square brackets:

F1 := -e*u*w-d*u+h*u+r1*z:

F2 := e*u*w-a*v-s*v*y+(r2+r3)*z:

F3 := k*v-u*w:

F4 := r3*z-b*y-(r2+r1)*y:

F5 := s*v*y-r*z:

solve({F1, F2, F3, F4, F5}, {u, v, w, y, z});

restart;

with(plots):

pr := 0.72: p := 0: n := [2, 3, 4, 5]: s := 1: a := 0.2: b := 1: L:=[red,blue,green,gold]:

for j to nops(n) do R1 := 2*n[j]/(1+n[j]); R2 := 2*p/(1+n[j]);

sol1 := dsolve([diff(diff(diff(f(eta), eta), eta), eta)+f(eta)*(diff(diff(f(eta), eta), eta))+R1*(1-(diff(f(eta), eta))^2) = 0, diff(diff(theta(eta), eta), eta)+pr*s^f(eta)*(diff(theta(eta), eta))+R2*pr*s*(diff(f(eta), eta))*theta(eta)+2*(a*(diff(f(eta), eta))+b*theta(eta))/(n[j]-1) = 0, f(0) = 0, (D(f))(0) = 1+b*((D@@2)(f))(0), (D(f))(5) = 0, theta(0) = 1+s*(D(theta))(0), theta(5) = 0], numeric, method = bvp); fplt[j] := plots[odeplot](sol1, [eta, diff(diff(f(eta), eta), eta)], color=L[j], axes = boxed);

tplt[j] := plots[odeplot](sol1, [[eta, theta(eta)]], color=L[j], axes = boxed) end do:

plots:-display([seq(fplt[j], j = 1 .. nops(n))]);

plots:-display([seq(tplt[j], j = 1 .. nops(n))]);

Example:

P:=2+x^5+x*y^3+x^2*y^4:

map(t->degree(t,x)!*degree(t,y)!*t, P);

                                              48*x^2*y^4+120*x^5+6*x*y^3+2

In this example, you can also use oblique quotes (this method is also suitable for earlier versions of Maple):

`4 + sqrt(9) = x - 1`:

parse(%);

%;

x=solve(%);

You can also delay the calculation of  sqrt(9)  by using of direct quotes:

4 + 'sqrt(9)' = x - 1;

%;

x=solve(%);

Your code works correctly when to end each line with a semicolon or a colon:

schro := {diff(psi(x), x, x)-(alpha*x^4+x^2-energy)*psi(x) = 0};

ic := {psi(3) = 0, (D(psi))(3) = 1};

schro1 := subs(energy = 3.30687, alpha = .1, schro);

soln1 := dsolve(schro1 union ic, {psi(x)}, type = numeric);

with(plots):

odeplot(soln1, [x, psi(x)], -3 .. 3);

Example (vector-row of plots):

plots[display](<plot(sin(x), x=-Pi..2*Pi, color=red) | plot(cos(x), x=-Pi..2*Pi, color=blue)>);