First, here's my entry for the basic question: 

SubsetPairs1a:= (A::list(set), B::list(set))->
local a, j, b;
    [for a in A do [for j,b in B do if a subset b then j fi od] od]
:   

(This procedure is threadsafe.) On repeated testing with high iterations, this has been about 15% faster for me than vv's entry. I haven't yet fully analyzed Tom's entry.

This information can be easily displayed as a bipartitite graph. Here's a procedure that does essentially the same thing as the above except that it formats the information in a way that's easier to use for the display purpose:

SubsetPairs1b:= (A::list(set), B::list(set))->
local i, a, j, b;
    {for i,a in A do 
        [for j,b in B do if a subset b then {-i,j} fi od][] 
     od
    }
: 

To get the plot, do

E:= SubsetPairs1b(ans6, ans7);
GT:= GraphTheory:
plots:-display(
    GT:-DrawGraph(
        GT:-Graph(E), layout= bipartite,
        stylesheet= [
            vertexcolor= yellow, 
            vertexfont= [TIMES,BOLD,18], 
            vertexpadding= 2
        ]
    ), 
    size= [1200$2]
);

There are many options available to make this plot look better, such as color and a more-readable font for the vertex labels.

Is there any reason that you don't want to use 

is(x^2 + y^2 >= 0) assuming x::real, y::real

?

The command is is the primary command for assessing the truth of relations containing assumed variables. However, if you prefer, you can make verify work with any boolean evaluator, including is. If you want this, let me know. 

To obtain a real solution, assign the results of solve to S (or any otherwise unused variable). Then do

(FreeV, BoundV):= selectremove(evalb, S):
FreeV:= lhs~(FreeV) =~ 0:
eval(BoundV, FreeV) union FreeV;

The 0 could be replaced by any number or any set or list of 7 numbers.
 

I believe that you're using Maple 2015 (because you've said so many times). You should always use the check boxes to put your Maple version in Question headers.

In Maple 2015, the precedence of the inert % operators is set incorrectly. This has been corrected in Maple 2020.

After you put M=0.5 in params, you get the error message

initial Newton iteration not converging

What I'm about the tell you applies only to this specific error message (and only if it starts with "initial"^[*1]). The technique is called a continuation parameter. You include a coefficient (which I'll call C here) in certain terms in the ODE system. You often need to experiment to find the right terms to use. C will vary from 0 to 1 during the initial phase of dsolve's solution process; so when C=0 it's as if the specified terms aren't there, and when C=1 ​​​​​​it's as if C wasn't there. Thus, it makes sense to include C in the complicated nonlinear terms. 

For this specific problem, I put C in every nonlinear term. Then I included option continuation= C in the dsolve command. It worked.

I also took out maxmesh=5000. Why would you impose a resource limitation (that may not even be needed) when you haven't even solved the problem yet?

[*1] Thus you need to be careful when discussing these errors here to distinguish between "Newton iteration not converging" and "initial Newton iteration not converging". Both of these occur quite often when solving boundary-layer BVPs.

The index (the part in square brackets) of a function call is accessible inside the function via the special syntax op(procname). So this works:

div5:= proc(f) local x; unapply(f(x)/5, x) end proc:
f:= x-> x^2+op(1, procname):
div5(f[5]);

The 1 is to select the first index, in case there are more than one. If there will always be exactly one, you can use simply op(procname).
 

My, what a trivial procedure, hardly worth naming, let alone writing a help page for. Anyway, here it is:

Trapz:= (X, Y)-> (X[2..]-X[..-2]).(Y[2..]+Y[..-2])/2:
trapz:= proc()
local X,Y;
    if nargs=1 then Y:= args[1]; X:= Vector([upperbound(Y)][1], i-> i)
    else X:= args[1]; Y:= args[2] 
    fi;
    if (local d:= rtable_num_dims(Y)) = 1 then Trapz(X,Y) 
    elif d=2 then Vector[row]([upperbound(Y)][2], i-> Trapz(X, Y[..,i]))
    fi
end proc:

I didn't handle the case where Y has more than 2 dimensions.

You can enter the system, numerically solve it, and plot the solution in 3D like this:

sys:= { 
    diff(x(t),t) = cos(x(t)), x(0) = 1, 
    diff(y(t),t) = sin(z(t)), y(0) = -1,
    diff(z(t),t) = z(t)*y(t), z(0) = 2
}:
sol:= dsolve(sys, numeric):
plots:-odeplot(sol, [x,y,z](t), t= 0..5);

If you want to solve it with Picard iteration, you'll need to code that yourself. Here's the plot produced by the above:

Do this:

f:= (x-1)/(x+2);
s:= convert(f, FormalPowerSeries);
c:= unapply((S-> (op(1,S),[op([2,1],S)]))(indets(s, specfunc(Sum))[]));
L:= limit(abs(c(k+1)/c(k)), k= infinity);
solve(L < 1, x); 
                   (-2,2)

Notes:

1. None of your simplify commands were needed, but neither would it be wrong for you to put them back.
2. Your line sum(k-> 1/k^2, k= 1..1000) is nonsense, and its result, 1000*f, is also nonsense. The fact that this doesn't give an error message should be considered a bug in Maple.
3. The line the turns the general term of the series into the function c is now fully automatic, no need to copy-and-paste.
4. The expression for the ratio test is abs(c(k+1)/c(k)). The extra x powers that you used are mathematical errors.
5. The final result is given in interval notation, the open interval from -2 to 2.
 

Here is how you can "automatically" add a plot to an existing plot:

AddPlot:= (p, P)-> plots:-display(P, p, _rest):
f:= (x-1)/(x+2):
p1:= plot(f, x= -3..3, discont, view= [-3..3, -5..5]):
p2:= AddPlot(
    plottools:-circle([0,0], 1, color= blue),
    p1,
    scaling= constrained
);

The following simple code is probably adequate for your needs. It could be made slightly more efficient by pre-sieving for the primes.

MyProperty:= proc(q::prime, P::name) 
local p:= ithprime(q - numtheory:-pi(q));
    if p >= q then P:= p; true else false fi
end proc
:    
MySeq:= proc(n::posint)
local k:= 0, R:= Vector(n), p, q:= 1;
    while k < n do 
        q:= nextprime(q);
        if MyProperty(q, 'p') then
            k:= k+1;
            R[k]:= [p,q];
        fi
    od;
    [seq](R)
end proc
:
S:= MySeq(9);
S := [[2, 2], [13, 11], [17, 13], [29, 17], [31, 19], [43, 23], 
  [67, 29], [71, 31], [97, 37], [107, 41]]

(p_seq, q_seq):= map(op~, [1,2], S)[];
   p_seq, q_seq := [2, 13, 17, 29, 31, 43, 67, 71, 97, 107], 
     [2, 11, 13, 17, 19, 23, 29, 31, 37, 41]

In Maple, type checking can only be done at run time. This is mostly because types can be arbitrarily complex, and checking them may involve execution of an arbitrary amount of Maple code. Of course, for the example that you show, it'd be theoretically possible before runtime; but I'm sure that you realize that your example is extremely superficial. 

The code below allows you to set the values of the parameters with sliders and have the plot update automatically. For some values of the parameters, there will be warning messages such as "cannot evaluate the solution further (left/right) of ... probably a singularity." These warnings do not prevent you from proceding with the Explore. If you want, the field arrows could be added.

restart:.
A:= a_1*c^3 - b_1:  B:= 3*a_1*c^2 - 3*b_1/c:  C:= d*c - e:
sys:= {diff(x(t), t$2) = (B*x(t)^2+C*x(t)-e)/A, x(0)=f, D(x)(0)=g}:
P:= [a_1, c, b_1, d, e, f, g]:
sol:= dsolve(sys, numeric, parameters= P):
MyPlot:= proc(P)
    sol(parameters= P);
    plots:-odeplot(
        sol,
        [x(t),diff(x(t),t)],
        t= -3..3,
        view= [-6..6, -6..6],
        thickness= 2, axes= boxed, color= red, scaling= constrained
    )
end proc
:
Explore(MyPlot(P), parameters= (P=~ -5..5));

subsindets(B, string, parse)

As Tom admits, explicitly entering the expression to be integrated is "ugly", and I find Axel's modification to be only a slight improvement. And that ugliness is made worse by the fact that the OP actually has a large number of such sets of integrals to do that I removed from the posted worksheet because the syntax issue is identical even though the integrands are not.

I like to build the expression to be integrated directly into the ODE system, like this:

restart;
eq1:=diff(f(y), y$4)+Uhs*diff(E(y),y$3)-(diff(f(y), y$2))+(diff(theta(y), y$1))= 0:
eq2:=diff(theta(y), y$2)+(diff(f(y), y$2)+1)^2+1+diff(theta(y),y$2) = 0:
E:=y->zeta*(cosh(k/2*(h1+h2-2*y)))/(cosh(k/2*(h1-h2))):
bcs:=f(h1) = -(1/2)*(Q-1-d), f(h2) = (1/2)*(Q-1-d), (D(f))(h1) = -1, (D(f))(h2) = -1,theta(h1) = 0, theta(h2) = 1:
epsilon1:=0.1:d:=1:omega:=Pi/6:h1:=-(1+epsilon1*sin(2*Pi*x)):h2:=d+epsilon2*sin(2*Pi*x+omega):F:= Q-1-d:epsilon2:=0.5:x:=1:
alpha:=Pi/6:
de:=eq1,eq2,bcs:
d1 := subs(Uhs =-2, zeta=3,k=1,{de}):
param:= {Uhs =-2, zeta=3,k=1}:
P1:= eval(diff(f(y), y$3)+Uhs*diff(E(y),y$2)-(diff(f(y), y$1)+1)+(theta(y))+sin(alpha),param):
ec:=0.5:
for Q from -3 to 3 by ec do 
    F2[Q]:= dsolve(
        d1 union {diff(J(y),y)=P1, J(h1)=0}, 
        numeric, maxmesh= 25500, continuation= lambda1
    ) 
end do:
for Q from -3 to 3 by ec do 
    P3[Q]:= eval(J(y), F2[Q](1)) - eval(J(y), F2[Q](0))  
end do;
                   P3[-3] := 2.30876985489011
                  P3[-2.5] := 1.36946865887356
                 P3[-2.0] := 0.460748411800256
                 P3[-1.5] := -0.417186862674401
                 P3[-1.0] := -1.26412891537822
                 P3[-0.5] := -2.07986512872413
                  P3[0.] := -2.86417840284371
                  P3[0.5] := -3.61684700714295
                  P3[1.0] := -4.33764445144679
                  P3[1.5] := -5.02633933429334
                  P3[2.0] := -5.68269517762163
                  P3[2.5] := -6.30647027769476
                  P3[3.0] := -6.89741755441151