The following technique will work regardless of the length of your code and regardless of whether it's in a Code Edit Region or a 1D-input (aka Maple Input) execution group.

# Here is some code in which I've intentionally included an error. Note that the first and last lines are 
# commented out.

# S:= "

P:= proc()
local j;
    for j to 9 do
        x:= x % 1
    od
end proc:

# ": 
Error, missing operator or `;`

# We got one of those typical error messages for which it's nearly impossible to locate its place in the code
# if the code were lengthy.
#
# Now I remove the # from the first and last lines. Thus, the entire code is put into a string. This is done in-
# place in practice; no need to copy-and-paste.

S:= "

P:= proc()
local j;
    for j to 9 do
        x:= x % 1
    od
end proc:

":

# The parse command will give more information than the previous error message:
parse(S, statement);
Error, incorrect syntax in parse: missing operator or `;` (near 55th character of parsed string)

# By using the information provided by parse, the following utility will help pinpoint the error.

HighlightError:= proc(S::string)
local m, p, nl1, nl2;
    try parse(S, ':-statement')
    catch "incorrect syntax":
        (m,p):= lastexception[[-2,-1]];
        nl2:= SearchText("\n", S, p..-1);
        nl2:= `if`(nl2=0, -1, p+nl2-1);
        nl1:= StringTools:-SearchAll("\n", S[..p]);
        nl1:= `if`(nl1=(), 0, [nl1][-1]);
        printf("%s%s^ ######%s", S[..nl2], cat(" "$(p-nl1-1)), S[nl2..]);
        error
    end try
end proc
:
HighlightError(S);


P:= proc()
local j;
    for j to 9 do
        x:= x % 1
              ^ ######
    od
end proc:


Error, (in HighlightError) incorrect syntax in parse: missing operator or `;` 
(near 55th character of parsed string)

Good Question. Vote up.

[Some of the points that I make below were also made by vv and acer while I was typing this. It was all written before reading those Answers, so please read it in that context.]

As I guess you've noticed, I is quite a strange object in Maple. If a Maple programmer wanted to create an object similar to I, it would be impossible to do using purely Maple code, although an alias would come close. This is different than say, Pi: If you wanted to make another constant similar to Pi, it'd be easy to do.

#First let's compare I with a "normal" constant like Pi:
OpTable:= (Op, R, C)-> 
    DataFrame(rtable(map(Op~, R, C)), ':-rows'= R, ':-columns'= C)
:
OpTable(type, [I, Pi], [symbol, atomic, complexcons]);
                      [    symbol  atomic  complexcons]
                      [                               ]
                      [I   false    true      true    ]
                      [                               ]
                      [Pi   true    true      true    ]

# I doesn't appear in the predefined global sequence constants:
constants;
               false, gamma, infinity, true, Catalan, FAIL, Pi

# Although the below doesn't answer the question of design (i.e., Why was it done that 
# way?), any such question of implementation can ultimately be answered by
# disassemble, addressof, and kernelopts(dagtag):
I_addr:= disassemble(addressof(I));
                      I_addr := 6, 18446744073709551617

kernelopts(dagtag= I_addr[1]), disassemble(I_addr[2]);
                                COMPLEX, 2, 1

kernelopts(dagtag= 2);
                                   INTPOS

addr:= disassemble(addressof(3+5*I));
            addr := 6, 18446744073709551621, 18446744073709551625

[disassemble]~([%[2..]]);
                              [[2, 3], [2, 5]]

In conclusion, I is like an alias: It's directly replaced by Complex(0, 1) without requiring evaluation. I'd guess that it's only done this way for historical and accidental reasons, and that it'd be done differently if Maple were being rewritten today from scratch.

Note that Complex(a, b), where a and b are manifest real constants, is not a function (in the ordinary sense that that word is used in Maple). Rather, it's one of the fundamental and directly-hardware-represented data structures, as shown by the disassemble and dagtag above. As shown above, that data structure has separate pointers to a and b except when a is 0.

The GUI (display) and kernel ("mathematical engine", which stores the variable values) are almost totally independent processes. There is no situation in Maple (nor has there ever been) where you can move to a previously executed portion of a worksheet and have the variables return to the values that they had when that portion was originally executed. So, this has nothing to do with Explore; it's just the way the worksheet interface has always worked.

One thing you can do is set the desired variable values in the same execution group that contains the Explore.

You can't use D as a parameter because it's the name of the differential operator used for the boundary conditions. This is the proximate cause of your error. You can change parameter D to D1 or whatever.

After that, there is a problem that a=0 reduces the differential order of the system to 3, so for this case only you need to remove one boundary condition.

After that, there is a problem that your 1st ODE needs to be explicitly solved for its highest-order derivative, diff(F(eta), eta$2). You need to do this either by hand (trivial), or with solve. Then you need to choose either the positive or negative branch of the square root.

After that, that first ODE is singular at eta=0 because the solution from the last paragraph has eta in the denominator. You can get around this by using option method= bvp[midrich] or method= bvp[middefer] in the dsolve command.

After that, there is a problem with negative values under the square root. There's no easy fix for this. Perhaps your parameter values are incorrect and you need to choose some such that that doesn't happen.

Keywords must be global, so they must be passed internally as ':-ic'= ic to avoid the name conflict with the parameter of the same name.

I wrote a threaded version of SubsetPairsSparse, and I achieved significant time savings in some cases.

restart:
SubsetPairsSparseThreaded:= proc(A::list(set), B::list(set))
local
    j,
    T:= table([for local x in `union`(A[],B[]) do x= table() od]),
    doB:= j-> (for local x in B[j] do T[x][j]:= () od),
    doA:= a-> `intersect`(for local x in a do T[x] od)  
;
    Threads:-Seq(doB(j), j= 1..nops(B));
    T:= {indices}~(T, 'nolist');
    Threads:-Map(doA, A)
end proc
:
SubsetPairsSparse:= proc(A::list(set), B::list(set))
local    
    a, b, x, j,
    T:= table([for x in `union`(A[],B[]) do x= table() od])
;
    for j,b in B do for x in b do T[x][j]:= () od od;
    T:= {indices}~(T, 'nolist');
    [for a in A do `intersect`(for x in a do T[x] od) od]
end proc
:
#My test case:
S:= {seq}(seq(a[i,j], i= 1..5), j= 1..5):
A:= [seq](combinat:-randcomb(S, rand(1..25)()), k= 1..2^13):
B:= [seq](combinat:-randcomb(S, rand(1..25)()), k= 1..2^13):

#Size of test case:
ProbSize:= (A::list(set), B::list(set))->
    [nops~, (add@nops~)~, (nops@`union`@op)~]([A,B])
:
ProbSize(A,B);
           [[8192, 8192], [106237, 106722], [25, 25]]

#Time both and compare results:
L1:= CodeTools:-Usage(SubsetPairsSparseThreaded(A,B)):
memory used=1.03GiB, alloc change=16.00MiB, 
cpu time=9.31s, real time=2.73s, gc time=734.38ms

L2:= CodeTools:-Usage(SubsetPairsSparse(A,B)):
memory used=1.01GiB, alloc change=-8.00MiB, 
cpu time=4.89s, real time=4.42s, gc time=1.09s

evalb(L1=L2);
                              true

kernelopts(numcpus);
                               4

As for the garbage collection settings, please test at least these two:

kernelopts(gcmaxthreads= iquo(kernelopts(numcpus), 2)):

kernelopts(gcmaxthreads= kernelopts(numcpus)):

I do not know how the competition for processor time between the garbage collector and Threads works. Hopefully this will give us some information about that.

Download PowSeriesTrigId.mw

Sorry, but MaplePrimes refuses to display inline any of my worksheets, no matter how simple. But you can download the above. Here's the text portion:

#Using power series to prove trigonometric identities

#Starting with the well-known Maclaurin series for sin(x) and cos(x), we prove sin(2*x) = 2*sin(x)*cos(x) with the help of Maple's symbolic summation.
(S,C):= convert~([sin, cos](x), FormalPowerSeries)[];

#The series' general terms can be rewritten like this:
(S,C):= Sum~((x^2)^k*(-1)^k/(2*k)!*~[x/(2*k+1), 1], k= 0..infinity)[];

#Multiply the series:
sin(x)*cos(x) = Sum(%Sum((-1)^j/(2*j)!*x/(2*j+1)*(-1)^(k-j)/(2*(k-j))!, j= 0..k)*(x^2)^k, k= 0..infinity);

#A closed form can be found for the inner sum, which is why I distinguished it as %Sum.
eval(%, %Sum= sum);
simplify(convert(%, factorial));

#Now compare that with the target:
sin(2*x)/2 = simplify(convert(sin(2*x)/2, FormalPowerSeries));

The type of an uninitialized module-wide variable with a type declaration is '`::`'(`local`, type) regardless of whether it's a local or export. For example

restart:
foo:= module() export n::integer; end module:
type(foo:-n, '`::`'(`local`, type));                
                              true

So, if you always give module-wide variables a type declaration (which is a good idea), this method can be used to check for them. Since anything and nothing are valid types, this can always be done. A variable declared nothing can never be assigned a value---it must remain purely symbolic---under the strict type checking imposed by kernelopts(assertlevel= 2).

Maple's type boolean is a larger class than you think; it's not simply identical(true, false). See ?type,boolean. There is a type truefalse, wihich is precisely identical(true, false). But Maple's true and false are predefined symbolic (i.e., nonnumeric) constants. In C#, are they numeric or symbolic constants? If numeric, then Is your intention to create a C# program that'll be called from Maple? Then you won't be able to directly pass Maple's true or false. But you could use a wrapper such as

TF:= (tf::truefalse)-> `if`(tf, 1, 0):

The wrapper would be used in Maple to pass such a value to the C# program.

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