Use  plottools[getdata]  and  plots[display]  with  insequence  option.

Example:

A:=plots[implicitplot](x^5+y^5-x*y=1, x=-3..3,y=-3..3, gridrefine=2):
B:=plottools[getdata](A):
C:=convert(B[3],listlist):
plots[display](seq(plot(C[1..n]), n=1..nops(C), 3), insequence);

You can use  seq  command:

DataName:=Array(["Data1", "Data2", "Data3"]):
DataSet:=[seq(i=DataName[i], i=1..numelems(DataName))];

                      DataSet := [1 = "Data1", 2 = "Data2", 3 = "Data3"]

Edited.

In this line 

fplt[j] := plots[odeplot](sol2, sol1, [[eta, diff(f(eta), eta)], [eta, diff(f(eta), eta)]], color = K1[j], axes = boxed);

should be one thing  only  sol2  or  sol1

Combine these graphs on one plot you can by  plots[display]  command.

By this way you can build angles on the plane anywhere and any values:

restart;
with(plottools): with(plots):
A:=line([0.5,0.5],[1.5,0.5], color=red, thickness= 3):
B:=rotate(A, 9*Pi/4, [0.5,0.5]):
C:=plot([0.5+(0.1+0.015*t)*cos(t),0.5+(0.1+0.015*t)*sin(t), t=0..9*Pi/4], color=black):
P:=curve([[0.655,0.6],[0.647,0.647],[0.68,0.625]]):
T:=textplot([0.45,0.45, 9*Pi/4], color=blue):
display(A, B, C, P, T, scaling=constrained, view=[0..1.7,0..1.5]);

Edited. Arrow added at the end of the arc, which shows the direction of the rotation.

If you want to immediately build a few graphs for specific values of the parameter m, and of specific colors, you can do so:

M:=[0.5, 2, 5, 8, 10]:  Color:=[red,green,blue,gold,black]:
plot([seq(-arctan(2*m*x/(1-x^2)), m=M)], x=0..10, color=Color,  numpoints=1000, discont, scaling=constrained);

All graphs have finite discontinuities at  x = 1  ( limit(f(x), x=1, left)=-Pi/2  and  limit(f(x), x=1,right)=Pi/2 )

Carl, you have incorrectly entered the expression for the function.

Pretty good approximation gives the model function with  arctan :

x_val:=<250,300,350,397,451,497,547,593,647,691,745,788,840,897>:

y_val:=<0,0.5,2,6.3,23.2,48.7,71.2,83.4,90.1,92.8,94.7,95.7,96.9,97.8>:

A:=plot(x_val, y_val, style=point, symbol=solidcircle, color=blue, symbolsize=12):  # Plot of initial data

f:=Statistics[NonlinearFit](a+b*arctan(c*x+d),x_val, y_val, x, initialvalues=[a = 50, b = 33, c = 0.02, d=-12]);  # Model function

B:=plot(f, x=x_val[1]..x_val[-1], color=red): # Plot of the model function

plots[display](A, B);

Addition: If the values of the y-variable are the percentages, then we can simply scale these values. The graph will have the same form:

x_val:=<250,300,350,397,451,497,547,593,647,691,745,788,840,897>:
y_val:=1/100*<0,0.5,2,6.3,23.2,48.7,71.2,83.4,90.1,92.8,94.7,95.7,96.9,97.8>:
A:=plot(x_val, y_val, style=point, symbol=solidcircle, color=blue, symbolsize=12):  # Plot of initial data
f:=Statistics[NonlinearFit](a+b*arctan(c*x+d),x_val, y_val, x, initialvalues=[a = 50/100, b = 33/100, c = 0.02, d=-12]);  # Model function​​​​​​​
B:=plot(f, x=x_val[1]..x_val[-1], color=red): # Plot of the model function
plots[display](A, B);

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](%);