The command lcm is primarily intended for polynomials, not integers. So lcm(3,k) is 3*k since k is a polynomial. Then the sum becomes sum(3*k, k= 0..3). The command for integers is ilcm, for which ilcm(3,k) will remain unevaluated until k has a numeric value.

Another solution is to use add instead of sum so that lcm(3,k) isn't evaluated until k has a numeric value. You should usually use add instead of sum for sums of a specific number of terms.

Any two plots with the same number of dimensions, regardless of package, can be merged with plots:-display.

plots:-display([plotA, plotB]);

If you use floating-point numbers in M by changing 2 to 2., then the unapply way is three times faster.

restart:
st:= time():
y := int(1/(-0.4016e-1*m^(2/3)-0.211e-3*m^(5/3)), m):
f:= unapply(abs(y), m):
n := 500: ## sample size
M := <seq(2.*idx/n,idx=1..n)>: ## m
Y := f~(M)+~Statistics:-Sample(Normal(0,3), n)^+: ## signal + noise
time()-st;

restart:
st:= time():
y := int(1/(-0.4016e-1*m^(2/3)-0.211e-3*m^(5/3)), m):
f:= (x) -> abs(subs(m=x,y)):
n := 500: ## sample size
M := <seq(2.*idx/n,idx=1..n)>: ## m
Y := f~(M)+~Statistics:-Sample(Normal(0,3), n)^+: ## signal + noise
time()-st;

                             0.156
                             0.547

However, I can't yet explain this phenomenon.

rsolve({x(n) = sqrt(7+x(n-1))/x(n-2)/9, x(1) = 2.0, x(0) = .25}, {x(n)}, makeproc)(852);

     0.1675558443

Another way: This is the most basic approach: a recursive procedure:

proc(n)
option remember;
     `if`(n=0, .25, `if`(n=1, 2., sqrt(7+thisproc(n-1))/thisproc(n-2)/9))
end proc
(852);

     0.1675558443

Another way: This is a little more efficient because it avoids checking the `if` conditions.

X:= proc(n) option remember; sqrt(7+thisproc(n-1))/thisproc(n-2)/9 end proc:
X(0):= .25:  X(1):= 2.0:
X(852);

Another way: This is the straightforward for-loop approach:

X[0]:= .25:  X[1]:= 2.0:
for n from 2 to 852 do X[n]:= sqrt(7+X[n-1])/X[n-2]/9 end do:
X[852];

Another way: This is more complicated to understand, but it is extremely efficient both in time and memory:

evalhf(
     proc(n, X0, X1)
     local x0:= X0, x1:= X1, x:= x1;
          to n-1 do
               x:= sqrt(7+x1)/x0/9;
               x0:= x1;
               x1:= x
          end do
     end proc
     (852, .25, 2.)
);

Another way: This is also extremely efficient in memory, and it's faster than the above if you're going to call it enough times to amortize the time taken for the initial compilation.

X:= proc(n::posint, X0::float[8], X1::float[8])
option autocompile;
local x0:= X0, x1:= X1, x:= x1;
     to n-1 do
          x:= sqrt(7+x1)/x0/9;
          x0:= x1;
          x1:= x
     end do
end proc:
X(852, .25, 2.);

This last procedure can very quickly calculate for n values in the millions.

You can't share equation labels between worksheets (as far as I know), but you can share variables, and anything that you can assign to a label could just as well be assigned to a variable. Go to the Tools -> Options -> General menu. The first item is "How should Maple handle the creation of new Math Engines?" Select "Ask each time a new document is created." (Actually, it will ask you each time any document, new or old, is opened.) Click "Apply to Session" or "Apply Globally," whichever suits you.

Now you're back in your worksheet. On the very top bar of the window, in the middle, it tells you the server number of that worksheet. Remember that number. Now open the worksheet that you want to share with. If it's already open, close it and reopen it. It will ask you to select a server, giving you a list of active servers. Choose the appropriate number.

Now every name, variable, matrix, etc. is known by the same name in both worksheets. There's no limit to the number of worksheets that can share a server.

In all cases, we have

(light ray vector) = (torus vector) - (light source vector).

Having the light source at the origin obscured the simplicity of that by makiing the light-ray vector equal to the torus vector. What matters (in all cases) is where the light-ray vector intersects the ellipsoid.

Here's the modified code. I removed the obfuscating factors like the title so that I could concentrate on the mathematics. You'll need to put those back.

restart:
macro(P= plots, LA= LinearAlgebra):
Ellipsoid:= V-> V[1]^2/32 + V[2]^2/18 + V[3]^2/12 = 1:
DispEllipsoid:= P:-implicitplot3d(
     Ellipsoid([x,y,z]), ([x,y]=~ -10..10)[], z= 1.25..5,
     style= surface, color= yellow, grid = [30$3]
):
TorusVorig:= <(1+.25*cos(v))*~[cos,sin](u)[], .6+.25*sin(v)>:
DispTorusVorig:= plot3d(TorusVorig, ([u,v]=~ 0..2*Pi)[], grid= [15$2]):
for X from -.8 by .1 to .81 do
     NormLightV:= LA:-Normalize(TorusVorig - <X,0,0>, Euclidean);
     Vadj:= solve(Ellipsoid(w*NormLightV), w);
     DispTorusVproj:= plot3d(
          Vadj[`if`(evalf(eval(Vadj[1]*NormLightV[3], [u,v]=~ Pi)) > 0, 1, 2)]*NormLightV,
         ([u,v]=~ 0..2*Pi)[]
     );
     Disp[X]:= P:-display(DispTorusVorig, DispEllipsoid, DispTorusVproj)
end do:
P:-display([entries(Disp, nolist, indexorder)], insequence, labels= [x,y,z], axes= boxed);

Perhaps you want a solution in the least-squares sense, i.e., a solution that minimizes the 2-norm of the residuals. That can be obtained like this:

Sol:= LinearAlgebra:-LeastSquares({seq(A[k], k= 1..20)}, {seq(C[k], k= 1..20)}, optimize);

Sol := {C[1] = 911044657906/553895179357, C[2] = 1720500130/3953256633, C[3] = 11884850/3953256633, C[4] = -50258975515/68523114972, C[5] = -38617059741145/13293484304568, C[6] = -7250038353509/6646742152284, C[7] = 13536241085573/9970113228426, C[8] = 18422595575075/39880452913704, C[9] = -21047503748893/9970113228426, C[10] = 201112232989883/39880452913704, C[11] = -3190968892022/1661685538071, C[12] = 189490779930919/39880452913704, C[13] = 109099913215741/39880452913704, C[14] = -6046079594847/2215580717428, C[15] = -42386684518409/6646742152284, C[16] = 26949111613257/4431161434856, C[17] = -59925038327045/6646742152284, C[18] = 431312338718351/19940226456852, C[19] = -206193995747093/4985056614213, C[20] = 149086421730665/3323371076142}

evalf(Sol);

{C[1] = 1.644796149, C[2] = .4352108375, C[3] = 0.3006344162e-2, C[4] = -.7334601694, C[5] = -2.904961473, C[6] = -1.090765700, C[7] = 1.357681781, C[8] = .4619454953, C[9] = -2.111059651, C[10] = 5.042877357, C[11] = -1.920320553, C[12] = 4.751470108, C[13] = 2.735673876, C[14] = -2.728891594, C[15] = -6.377061656, C[16] = 6.081726430, C[17] = -9.015700768, C[18] = 21.63026281, C[19] = -41.36241806, C[20] = 44.85999857}

Use an animation with respect to w:

plots:-animate(
     plot3d,
     [[(5+w*cos(v))*cos(u), (5+w*cos(v))*sin(u), w*sin(v)], u= 0..2*Pi, v= 0..2*Pi],
     w= 0..3
);

Use a parametric plot, like this:

y := int(1/(-0.4016e-1*m^(2/3)-0.211e-3*m^(5/3)), m):
plot([y, m, m= 0..2], labels= ['y', m]);

Use plots:-spacecurve to plot the 2D curve with a z-coordinate of 0.

plots:-display(
     plot3d(exp(x*y), x= 0..1, y= 0..1),
     plots:-spacecurve([x, exp(x), 0], x= 0..1, color= red, thickness= 2)
);

A complete plot of exp(x+y+z) would require four dimensions. Using three dimensions, you'll have to settle for the function's level surfaces, which can be plotted via

plot3d(-x-y+ln(C), ...)

where C is the level value.

For a < b, Heaviside(x-a) - Heaviside(x-b) is 1 for a <= x < b and 0 otherwise (given your redefinition of Heaviside(0)). Therefore, the integral

Int(f(x), x= a..b)

is equivalent to

Int(f(x)*(Heaviside(x-a) - Heaviside(x-b)), x= -infinity..infinity).

Avoiding automatic simplifications is a very tricky thing, and I don't think that there's any perfect solution. Your case is especially tricky because you also want to preserve superfluous parentheses. I think that the best solution is offered by the InertForm package.

macro(IF= InertForm):

Enter the expression as below. The purpose of the ``s is to preserve the parentheses.

ex:=
     ``(IF:-Parse("F^5*alpha[5]+F^4*alpha[4]+F^3*alpha[3]+F^2*alpha[2]+F*alpha[1]")) +
     ``(IF:-Parse(
            "F^4*G*gamma[1]+F*G^4*gamma[3]+F^2*G*gamma[2]+"
            "F*G^2*gamma[4]+F*G*gamma[5]"
     )) +
     ``(IF:-Parse("G^5*beta[5]+G^4*beta[4]+G^3*beta[3]+G^2*beta[2]+G*beta[1]")) +  
     delta[0]
:

To view the expression, use

IF:-Typeset(ex, inert= false);

To use the expression mathematically, use

expand(value(ex));

It's not the complexity that matters. The problem is that (as far as I can tell) geometry:-draw can only plot within the window [-10..10, -10..10]. Your first line doesn't intersect that window.

As a workaround, use implicitplot:

with(geometry):
point(P1, [47+(38+22/60)/60, -(122+(43+4/60)/60)]):
point(P2, coordinates(P1) +~ [cos,sin](30*Pi/180)):
#Specify axes names as 3rd argument to `line` to avoid super-annoying dialogs!
line(L1, [P1,P2], [x,y]):
plots:-implicitplot(Equation(L1), x= -100..100, y= -100..100);

Update: The second sentence above is incorrect: The default view is [-10..10, -10..10], but that can be changed or nullified.

Change the line D:= A to D:= copy(A). The line D:= A doesn't mean what you think it does. It means that D will be A and will stay A until it is reassigned. Thus, any changes to D also change A.

There's a problem with how you entered T: You can't use square brackets for algebraic groupings. Enter T like this:

T:= (r,theta)-> 0.5*alpha+Sum(gamma*n*(r^n+R^(2*n)+r^(-n))*sin(n*theta), n= 1..infinity):

Use PDEtools:-dchange to make the change of variables:

PDEtools:-dchange({theta= arctan(y,x), r= sqrt(x^2+y^2)}, Diff(T(r,theta), y), [x,y]):
value(%);

Notice that Diff is with a capital D. This is necessary to delay the evaluation of the derivative until after the change of variables is applied. The value then performs that evaluation, which takes about three minutes, unfortunately.