Yes, I think that the defaults should be changed. Specifically, 

1. It's often difficult to appreciate the 3-dimensional nature of a plot without gridlines, though, unlike your example, I wouldn't necessarily include them on both the inner and outer surfaces unless they were differently colored (which IS fairly easy to do).
2. With the current defaults, the colors of the inner and outer surfaces are too similar to each other.

@rlopez Dr Lopez's suggestion can be given the full sign, primpart, sort treatment to put the results into your requested format, doing them all at once with this single command:

`[]`~(
    L,
    (sign*primpart@sort@lhs@Student:-VectorCalculus:-TangentPlane)~(
        (x - 1)^2 + (y - 3)^2 + (z - 5)^2 = 13^2, `<,>`~(L)
    ) =~ 0
);

I am suprised that this Student command doesn't check whether the specified point is actually on the specified surface. It'll happily return results for points that are not. (This is mentioned on its help page.)

@NIMA112 Why do you suggest this obscure numerical method that'll require much work to get coded into Maple? I'm not saying that it's not a good method for this problem---that I don't know. Do you have any reason to suspect that this method will work better than any of the numerous more-standard methods that are already coded?

@Art Kalb That's a good question, but I don't know the answer. It is perhaps a shortcoming or bug of expand.

@Laurenso Here is the version of the procedure with the on-the-sphere check included:

tangent_plane:= p-> local P:= <p>, R:= <1,3,5> - P;
    [p, `if`(R.R = 13^2, (sign*primpart@sort)(R.(<x,y,z> - P)) = 0, FAIL)]
:

An equivalent to map(expand@expand, N^%H) is 

(expand~@@2)(N^*)

If you want some variables to represent real numbers, then you must use assume or assuming. If you want them all to represent reals, then why not just use the ordinary transpose?

@Laurenso I just realized that my procedure doesn't check whether the given point is actually on the sphere. If you require that check, I can easily add it.

@Laurenso The command to see the procedure is showstat(geom3d::onobjl).

@2cUniverse I think that you've found some relevant info. It'll take me a few days to digest it. Thanks.

Regarding your Reply before the one above: As I said before, that idea is already implemented in the presented code. Once any element x of order k is found, I immediately use it to generate another totient(k) - 1 elements of order k. What I'm missing is a fast way to generate the initial elements x needed to get all elements of order k. Starting with elements of order Carmichael(n) isn't enough.

So, do you want all pairs of perpendicular diameters with integer endpoints? 

It looks like you're trying to solve a system of 41 equations with only 8 unknowns.

@dharr I concur: vb has been changed from 0.5 to 5. Since vb is the numerator of the arguments of arccos, that change would have a major impact on the chance of there being a real solution (since arccos(x) is not real for real x such that abs(x) > 1).

Using the original vb = 0.5, I get this solution:

sol := {A = -1.261525434, B = 1.459322610, C = 0.03167646665,
    omega = 2.083929566, x1 = -1.965825318, x2 = -0.8500707662}

@dharr If NLPSolve detected complex numbers, then its error message would directly mention that. Since the error says "non-numeric", the input very likely still contains either symbolic non-decision variables or items with inappropriate indices (such as 2[1]).

I'm not saying that the input doesn't generate complex values; I don't know whether it does. I'm just saying that that can't be the main issue.

As you said, the OP needs to upload the worksheet.

@dharr The proximal cause of the error is that SparseDirectMKL is declared local to LinearAlgebra (for some unknown and likely unintentional reason), but the global SparseDirectMKL is being passed by the user's command. If this alone is the problem, it's easy to workaround; but there's another problem also:

kernelopts(opaquemodules= false):
:-SparseDirectMKL:= LinearAlgebra:-SparseDirectMKL:
kernelopts(opaquemodules= true):
sol:= LinearAlgebra:-LinearSolve(A, b, method= SparseDirectMKL);
Error, (in LinearAlgebra:-LinearSolve) data incompatible with method = LinearAlgebra:-SparseDirectMKL