You were very close to the solution

n := `2`:
p := ListTools:-Search(n, lhs~(S));
[n, rhs(S[p]):-mu, rhs(S[p]):-sigma]
                               3
             [2, 1208.001018040, 141.512246146314]

A solution which reuses the most part of what you did (this was my my constrraint)
 

restart:
with(ListTools):
C := x*u[]*u[1,2,2]^3*u[1,1,2];
indx_set := [ (`[]`@ op)~(indets(C, indexed))[] ];

deg      := degree~(C, [indets(C, indexed)[]]);
Indx_set := zip((u,v) -> u$v, indx_set, deg);

Flatten(convert(Indx_set,list));

ListTools:-Occurrences(1, %);
                        5

I'm not sure at all I understood correctly, but maybe this could help?
if I'm out of line just forget it
 

# If only the sums of terms with ranks 5, 10, 15, ...matter, what is the need to 
# compute U and V?

restart:
u := (a, b) -> exp(-(a-Lmu)^2)*exp(-(b-Wmu)^2);
llim := 1:
ulim := 20000:
step := 5:
nok := iquo(ulim-llim+1, 5);
        4000
W := CodeTools:-Usage( Vector(nok, i -> add(u(5*i, x__2)*u(x__4, 5*i))) );

memory used=19.89MiB, alloc change=36.00MiB, cpu time=313.00ms, real time=314.00ms, gc time=21.15ms

Download one_way.mw

This was intended to answer your last question, but it seems to have disappeared.
The procedure Circle in the attached file finds, if it does exist, the "tritangent" circle for arbitrary rational fraction P(x)/Q(x) where:

• P(x) is of degree at most 2
• Q(x) has degree 2 and two real roots

The procedure Circle solves the equation 

(x-p)^2+(f(x)-q)^2-r^2=0

whose numerator is a polynom of degree 6 in x.
The circle of centre [p, q] and radius r has 3 tangency point if this numerator is of the form 

(x-A[1])^2 * (x-A[2])^2 * (x-A[3])^2

After having equalized the coefficients of these two equations the procedure solves a system of 6 equations (p, q, r, plus  the abscissas A[1], A[2], A[3] of the three tangency points).


The file contains two "pathological" examples:

• In the first one, even if it seems visually that a solution exists, procedure Circle fails to find one.
  In procedure Circle change the line

  sol := {map(s -> if keep(s)  then s end if, sols)[]};

  by the line

  sol := sols;

  to see what happens: 

  • the function keep is aimed to discard solutions with nerative r, identical contact points or imaginary contact points

  • writting sol := sols instead will return the "unfiltered" solutions

• The second one is very interesting because it shows that for some functions f(x) with three branches, there can exist a circle which tangents a branch once and a second branch twice.
   

Another_way.mw

Kitonum's variant.
Use tubeplot instead of spacecurve:

• pros: doesn't require a fictitious shift of the ellipse (with is to be adjusted in an ad hoc way)
• cons: for scaling=unconstrained (the default value) the tube doesn't have a circular cross section but an elliptical one: if you increase the value of radius the ellipse will ressemble more to a short thin sylibder than to a line.

So make your choice

tubeplot([cos(t), sin(t), 1 - 12*cos(t) + 15*sin(t)], t = 0 .. 2*Pi, radius=0.05, color=red, orientation = [-15, 68, 5])

Seems there is no solution.

Solution without range constraints

N := 4:
sys := NULL:
for k to N while k <> j do 
sys := sys, 
       sinh(x[k]+I*a)^N*(product(sinh(1/2*(x[k]-x[j]-(2*I)*a)), j = 1 .. N))
       /
       (sinh(x[k]-I*a)^N*(product(sinh(1/2*(x[k]-x[j]+(2*I)*a)), j = 1 .. N))) 
end do:
sys := [sys]:

# check
# print~(sys):

fsolve(sys=~[-1$4], {x[1], x[2], x[3], x[4]});
                               
            {x[1] = -3.799588601 - 1.570796327 I, 

              x[2] = -3.799588601 - 1.570796327 I, 

              x[3] = 0.1196098914 - 1.570796327 I, 

              x[4] = 1.110549099 - 1.570796327 I}

# check solution (t[1]..t[4] are your original variables) 
# eval([seq(t[i], i=1..4)], %);

You cannot impose ranges the way you did.
The correct way is to write something like a1+b1*I .. a2+b2*I
For instance, with a1=-10 and a2=+10, b1=-1, b2=+1 one gets

fsolve(sys=~[-1$4], {x[1]=-10-1*I..10+1*I, x[2]=-10-1*I..10+1*I, x[3]=-10-1*I..+10+1*I, x[4]=-10-1*I..+10+1*I});
            {x[1] = 0.2343030315 - 0.9904555257 I, 

              x[2] = 0.9120742051 - 0.5349352395 I, 

              x[3] = 0.2227793496 + 0.3919592237 I, 

              x[4] = 0.2227793496 + 0.3919592237 I}

For the imaginary range you want there is unfortunately no solutions found (mayba you will find for other values of a1 and a2???)

fsolve(sys=~[-1$4], {x[1]=-10-I..10-0.5*I, x[2]=-10-I..10-0.5*I, x[3]=-10+0.5*I..+10+I, x[4]=-10+0.5*I..+10+I});

Remark
Depending on the way the system is written  fsolve will find different solutions.
Example:

sys_nd := map(simplify, numer~(sys)+~denom~(sys)):
fsolve(sys_nd);

            {x[1] = 0.7772186938 + 0.1314607556 I, 

              x[2] = 0.7772186938 + 0.1314607556 I, 

              x[3] = -0.7600984061 - 1.168722678 I, 

              x[4] = 0.7772186938 + 0.1314607556 I}

#check the solution
#eval([seq(t[i], i=1..4)], %);

But even formulated this way there is still nosolution found for the imaginary ranges you want:

fsolve(sys_nd, {x[1]=-10-I..10-0.5*I, x[2]=-10-I..10-0.5*I, x[3]=-10+0.5*I..+10+I, x[4]=-10+0.5*I..+10+I});

Shorter but less elegant than @MapleMathMatt 's, it requires a few knowledge on MathML syntax
 

sinx := `#mrow(mo("sin",mathcolor="red"),mo("(",mathcolor="red"),mo("x",mathcolor="red"),mo(")",mathcolor="red"))`:
plot( sin(x), x = 0 .. 2 * Pi, title = typeset( "Plot of %1", sinx ) );

If the space between "(" and "x") doesn't suit you, replace the definition of sinx by
 

sinx := `#mrow(mo("sin",mathcolor="red"),mo("(",mathcolor="red"),mo(" "),mo("x",mathcolor="red"),mo(")",mathcolor="red"))`:

A more formal way
 

restart;
with(codegen):

B := p -> x^4 - 4*x^3/(2 + p) - 6*(p - 1)*x^2/((3 + p)*(2 + p)) - 4*p*(p - 5)*x/((4 + p)*(3 + p)*(2 + p)):

# formal roots of B(p)

s := solve(B(p), x):
a  := [allvalues(s)]:

# The problem could end here, but plotting the solution a for different values of 
# p would require a lot of time for a has a very lengthy expression as cost 
# (codegen:-cost) demonstrates

cost(a[1]);
1246*additions + 6298*multiplications + 136*divisions + 162*functions

# the idea is to use the package codegen to rewrite a into a form which 
# minimizes the number of operations it contains.
#
# first step : transform each a[k] into a procedure with argument p

for k from 1 to 4 do
  F||k := makeproc(evalf(a[k]),p):
end do:

# second step : join the 4 procedures F1..F4

G := joinprocs([seq(F||k, k=1..4)]):

# third step : optimize the joint procedure G

H := optimize(G):

# cost of H (must be compared to 4 times the cost of a[1])

cost(H)
41 storage + 41 assignments + 107 multiplications + 90 additions + 15 divisions + 10 functions

# plot

CodeTools:-Usage(
  plots:-complexplot(
    [seq(H(q), q=0..5000)], 
    style = point, symbol=solidcircle, symbolsize=6
  )
);
memory used=170.05MiB, alloc change=32.00MiB, cpu time=1.97s, real time=1.80s, gc time=258.09ms

# for comparison

CodeTools:-Usage(
  plots:-complexplot(
    [seq(op(evalf(eval(a, p=q))), q=0..5000)], 
    style = point, symbol=solidcircle, symbolsize=6
  )
);
memory used=0.72GiB, alloc change=32.00MiB, cpu time=6.72s, real time=6.29s, gc time=624.36ms

As you can see in the attached file
exact_plus_codegen.mw
the formal approach is about 3 times faster than the numerical one and uses about one third of the memory this latter uses.
 

If you have written something like

a := 0.4:
:
:
:
plot(a, a=0..1)

and if you want ti unassign a before the plot command, just do

a := 0.4:
:
:
:
a := 'a':
plot(a, a=0..1)

This answer is based upon Maple 2015.2
If you have more recent version of Maple (which I don't at home) , DrawGraph proposes a lot feautures to customize the graph.
In any case, whatever the version, you can always modify the plot the way you want.

Here is an example.
Note that at some point I use GetVertexPosition to catch the position of the vertices: you can modify them as you want (SetVertexPosition).
As shown in tha attached file, you can also dilate a graph in a simple way by using  plottools:-transform.
 

small_vertex_names.mw

Put a DEBUG(); instrruction within your procedure Fig and execute it.
You will see that the source of the message is solve returning nothing.
To verify this put print(sol) after the solve command.

What is the purpose of this solve: to find the intersection between a couple of lines and an ellipse (at least for alpha =Pi/4).
In this case you can replace sol := solve({ell, OPQ}, {x, y}, explicit);  by this
 

__Z1 := [solve(OPQ, y)];
map(u -> eval(ell, y=u), __Z1);

# x coordinates of the 4 intersection points

sol_x := solve~(%);

# y coordinates of the 4 intersection points (if necessary)

sol_y := [ 
           map(u -> eval(__Z1[1], x=u), sol_x[1..2])[],
           map(u -> eval(__Z1[2], x=u), sol_x[3..4])[]
         ]:

# [x,y] coordinates of the 4 intersection points (if necessary)

sol := zip((u,v) -> [u,v], sol_x, sol_y)

The Statistics package has nothing to do this.
You will find in the attached file the solution for standard gaussian random variables (GRV) ALONE.
The method is classical and uses  the Cholesky decomposition of the correlation matrix.
A second exampleisprovided forGRVs with different variances.

If at least one variable is not GRV the notion of correlation is ambiguous, let me know if your are in this case

 

Download Correlated_gaussian_RV.mw

It's always possible to call Maple from any language which possesses a "call system" feature.
I don't have Fortran right now to verify this, but this should work:
 

call system ('\\...\maple file')   #Linux or Mac OS example

Remarks:

• On Windows platforms change maple into cmaple
• \\...\ is the absolute path to the maple executable file (it depends on your installation)
• file is the name ofthe mw file to run
• about maple/cmaple look here view.aspx

A few advices:

• remove all the plots unless you are fan of old-fashioned ascii-look plots
• remove warnings (maple's option -w)
• remove all interactive features in your worksheet
• as Maple will not communicate with your Fortran program, use a file mechanism to exchange informations:
  • write in a file named IN (for instance) the informations you want to pass to Maple
  • modify your Maple worksheet for it reads file IN and writes its results into a file named OUT whose format can be read by Fortran.
• some features are not supported by maple/cmaple : avoid using the Grid, Threads, DocumentTools and Maplets packages (probably not a comprehensive list)

I use to use this 'call system' strategy to call Maple from other languages (and reciprocally to call from Maple apps written in other languages, for Maple proposes a system/ssystem feature): a thing you have to be aware of is that this call is not instantaneous, so avoid repetitive calls and adapt your codes accordingly.

There are several issues to solve:

1. A level curve can be connected or not.
   Even if connected its drawing on some bounded domains may contains non connected pieces.
2. In the simplest case of a level curve made of a single piece in the drawing domain: do you want your data file to contain a set of points described in a a "natural order" (induced by the direction used to travel the curve) or doesn't it matter to you?
3. If there are multiple level curves and/or non connected level curves, I guess you would like to identify each branch one way or another. Thus your data file cannot only contain a N by 2 matrix: What structure did you think of?

Now a "Maple issue":

• contourplot builds a level curve made of a collection of segments [P, Q] where P and Q are two distinct points.
  Extracting those segments from the CURVES plot command shows they are not sorted in a "natural order" (see an illustration in the attached file).
  This is where point 3 above becomes important: "reordering" these segments seems extremely difficult (if possible?) in a general case.

In case the "order" matters, a proposition which doen't use contourplot is given in the attached file.

This file doesn't give a definite solution, just a few ideas.

ideas.mw

restart:
interface(version)
Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895

with(Statistics):
UseHardwareFloats :=  false:

# random construction of the kind of data matrix (named S) you have

C1 := RandomVariable(DiscreteUniform(1, 4)):
C2 := RandomVariable(DiscreteUniform(1, 9)):
R  := RandomVariable(DiscreteUniform(0, 100)):

N := 100:
S := < round~(Sample(C1, [N, 1])) | round~(Sample(C2, [N, 1])) | Sample(R, [N, 4]) >


# An example of a few combinations

Choice := [ [1, 3], [3, 4], [4, 6], [4, 8], [5, 10] ]:

# Extract all the couple of indices (cols 1 and 2) of S

cols12 := convert(S[.., 1..2], listlist):

# Search all the repetitions of each choice within the list cols12

occ := table([ seq(u = [ListTools:-SearchAll(u, cols12)], u in Choice) ])

# Extract the submatrix corresponding to each choice

table([ seq(u=S[occ[u], 3..-1], u in Choice) ])

# Or, extract the list of numbers corresponding to each choice

table([ seq(u=[entries(S[occ[u], 3..-1], nolist)], u in Choice) ])

PS: I used Maple 2015 and convert(.., listlist) is obsolete for newer version (sorry I do not have in mind the correct option to use to convert a matrix into a list of list, but type ?listlist to find it).

SubSample.mw