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.

What you call a perfect matching, I call a bijection. Here's a little procedure that finds all bijections between lists A and B:

Bijections:= proc(A::list, B::depends(And(list, satisfies(B-> nops(B)=nops(A)))))
     map[3](zip, `[]`, A, combinat:-permute(B))
end proc:

Bijections([1,2,3], [4,5,6]);

[[[1, 4], [2, 5], [3, 6]], [[1, 4], [2, 6], [3, 5]], [[1, 5], [2, 4], [3, 6]], [[1, 5], [2, 6], [3, 4]], [[1, 6], [2, 4], [3, 5]], [[1, 6], [2, 5], [3, 4]]]

If you have Maple 2016, then the following may be more efficient:

Bijections:= proc(A::list, B::depends(And(list, satisfies(B-> nops(B)=nops(A)))))
local b;
     seq(zip(`[]`, A, convert(b[],list)), b= Iterator:-Permute(B))
end proc:

Remove the line

xi:= k[1]*x-k[1]*c[1]*t ;  eta:= k[2]*x-k[2]*c[2]*t ;

Replace it with

PDEtools:-dchange({xi= k[1]*x-k[1]*c[1]*t, eta= k[2]*x-k[2]*c[2]*t}, Diff(u,t), [x,t]):
value(%);

In Maple, the end of a line is treated like a space, so there's no need for any line-continuation symbol. For human readability, the continued parts of the lines should be indented, but the code parser doesn't care. In Maple, your example command could be

A:= [
     3, 4, 5, 6, 6, 45, 37,
     5, 4, 67, 39, -967
];

In the extremely rare situation that you need to end a line at a place where a space wouldn't be allowed, end it with a backslash \. But I can't think of any situation where there's not a more-elegant way to do this.