Here are procedures for the primary and inverse functions:

Lehmer:= (S::And(list(nonnegint), satisfies(S-> max(S)=nops({S[]})-1)))-> 
   Rank(Iterator:-Permute([{S[]}[]]), S) - 1
:
LehmerInv:= (s::nonnegint, n::And(posint, satisfies(n-> s < n!)))->
   [seq(Unrank(Iterator:-Permute(n), s+1) -~ 1)]
:

Your fsolve is returning unevaluated for some loop iteration. This can mean several things, all of which commonly occur:

1. There is actually no solution.
2. fsolve thinks that there is no solution even though there actually is. In other words, some iterative process akin to Newton's method is diverging for every initial value that it tries.
3. fsolve knows that there is a solution, but it can't refine it to sufficient accuracy. In other words, the iterative process is converging, but there's some oscillation in the last few digits.

So, whenever fsolve is used programmatically, such as in a loop, you should check the form of the returned value (your S) to make sure that it isn't of the form fsolve(...) before that S is used in a further computation. Do something like this:

Failures:= table();
for i ... do
    ...; 
    S:= fsolve(...);
    if eval(S,1)::specfunc(fsolve) then 
       Failures[i]:= eval(S,1)
    else
       w(i):= eval(k, S)
    end if
end do: 
eval(Failures);

At the end, Failures will contain all the cases for which fsolve didn't converge. You may decide to analyze these further, or you may decide that the cases that did converge are sufficient for your purpose.

If you need further help with this, please post your complete code.

Here is a procedure to return the dependent variables from a PDE, ODE, or a set or list thereof:

depvars:= (pde::{algebraic, `=`(algebraic), {set,list}(algebraic, `=`(algebraic))})-> 
   indets(
      indets(convert(pde, diff), specfunc(diff)), 
      And(typefunc(name, name), Not(typefunc({mathfunc, identical(diff)})))
   )
:

And here is an implementation of is_solution_trivial:

is_solution_trivial:= (pde, sol::{`=`, set(`=`)})-> 
   evalb(eval(`if`(sol::set, sol, {sol}), depvars(pde)=~ 0) = {0=0}) 
:

I would guess that the following is a fundamental principle of engineering. I've never formally studied engineering, so I don't know whether it is. To me, it is just intuition (something that is obvious), having taking apart and rebuilt machines and other systems my entire conscious life. (I'd appreciate it if someone who knows would tell me whether it is indeed formally a fundamental principle of engineering.)

Let's say that you have two closely related systems, which I'll call W (for working) and N (for not working). They could be machines, programs, scientific theories, literary essays, biological systems, etc. So, W is working (i.e., it does what it was designed to do), but it needs to be modified so that it can do something else, or something additional. On the other hand, N has been designed and built based on W and it incorporates the needed changes, but it doesn't work.

My principle is 

1. Junk N. Discard it. It may be disassembled and small parts re-used. Otherwise, it should be forgotten, regardless of how much effort went into it.
2. Proceed by making small changes to W. Test after every small change to make sure that it's still working. Do not make any further changes until the current version works! This may require hundreds of small changes if you're working alone, and perhaps much more if you're working as a team.
3. (Optional step) Ocasionally remove from W parts that are no longer needed. The smaller size will make the work proceed faster.

Now, you may protest, "But a lot of effort and/or expense went into N!" That doesn't matter; it's more efficient and/or economical to work on W. Accountants refer to this as "sunk costs"; investment experts say "Don't throw good money after bad"; professors of creative writing say "Kill you darlings". If my principle were not true, then life on Earth would likely have started independently many times. Yet all evidence suggests that it hasn't.

So, in this thread, it seems clear that you have a W program and an N program. You've been trying to correct the N program. Stop that. Work on your W program.

{k = 1, n=1} satisfy your stated conditions. Other than that, I checked for n up to 54177 and found no other solution.

Like this:

convert(evalc(Im(D1C1)), exp)

I think that if you just view the indices by doing seq(nops(L)-i, i= 1..nops(L)), then you'll understand the reason that your second version doesn't work. If not, let me know, and I'll explain further.

By the way, you're finding the reverse of the list, not the inverse. The inverse is something else.

There's currently no predefined way to do that with the legend option. However, it wouldn't be difficult to simply overlay your plot with a legend (exactly as you show above) constructed with the plot and textplot commands. The overlaying is done with the display command. This requires knowing the approximate x and y coordinates of the placement of the legend, preferably the upper left and lower left corners. If you have foreknowledge of the full extent of the x and y axes, this is easy. If you don't have foreknowledge, it's possible to determine the axis extents of an existing plot programmatically.

Yes, it can be done with if ... then. The key is that all the code must be in a single execution group (in other words, a single command prompt). While typing code, to get a new line within the same execution group, use Shift-Enter instead of Enter (or carriage return). Example:

#There's no limit to the complexity of conditions that can be made with
#logical connectives such as `or`, `and`, etc.

if x > 0 then 
   y:= sin(x); #code line 1 to skip
   z:= cos(x)  #code line 2 to skip
   #... more code lines
fi; #or use "end if"

#Or, there can be alternative branches:

if x < 0 and a > 0 then
   y:= sin(x);
   z:= cos(x)
else
   y:= cos(x);
   z:= sin(x)
fi;

#If statements can be nested arbitrarily, and 
#"if ... then ... else if ... then ... fi fi" can 
#be contracted to "if ... then ... elif ... then ... fi".

The above are simplistic examples to get you started. There's actually easier ways to do each.

Here's a procedure that does it:

CheckProperty:= proc(n::And(odd, posint, Not(1), Not(prime)))
local nm1:= n-1;
   ormap(p-> irem(nm1, p[1]-1)=0, ifactors(n)[2])
end proc: 

Both numtheory:-factorset and NumberTheory:-PrimeFactors will do essentially the same job as ifactors here, but they do it by calling ifactors, which could be problematic for some very large  n (say n > 2^200).

The part that you want to extract is op([2,1], sol).

How about

plots:-textplot([0,0,typeset(Sigma)], font= [Times, 300], axes= none);

All that you need is

plot([seq([[k,0], [k, 3*k+2]], k= 1..6)]);

Making a plot rotate is quite simple. Like this (this code picks up immediately where your worksheet left off):

P:= %: #Save the static plot
Frames:= 20:
plots:-animate(
   plottools:-rotate, [P, theta], 
   theta= 0..2*Pi/5*(1-1/Frames), frames= Frames, paraminfo= false
);

The rotation doesn't need to be around the origin. It'd be trivial to rotate around any point.

The integrals have as their third argument the option AllSolutions, which latex can't handle. And what would you have it print for those anyway? So, you can workaround by doing

latex(subs(AllSolutions= (), sol))

The int in `latex/int` refers to the int and Int functions, not to "internal".