It can be done like this:

(n,k):= (6,2);
mul(add~(combinat:-choose([x||(1..n)], k))) = 0

The command plot3d is for plotting objects that are locally two-dimensional (what we commonly call surfaces) in three-dimensional space. What you want to plot in 3-D is an object that is locally one-dimensional (commonly called a curve). The command for that is plots:-spacecurve.

If you can figure out how to use it in this case, great, please post the results. If you have trouble figuring it out, I will provide further help.

I can't reproduce your situation; however, this may help, so please try it and let me know how it goes: A restart command should always be separated into its own execution group. The problems that can arise if you don't do this are often weird and difficult to explain and often involve the display of output.

Also: Although it is not an error and certainly not a source of any print problems, you should remove the line iter:= iter+1 from the loop. The for statement automatically increments its index variable for each iteration of the loop. (This is also true for the for statements in every other computer language that I know that has a for statement.) You should always strive to avoid making assignments to a for loop index inside the loop. It itself is not an error, but it is frequently a source of erroneous results.

Consider changing your style to something like this:

eps:= 1e-12; itermax:= 10; Error:= 10;  
for iter from 0 while Error>eps and iter<itermax do
   print("iteration", iter, "Error", Error); 
   Error:= Error*0.1;      
end do;

or this:

eps:= 1e-12; itermax:= 10; Error:= 10;  
for iter from 0 to itermax-1 while Error>eps do
   print("iteration", iter, "Error", Error); 
   Error:= Error*0.1;      
end do;

These are more-normal styles than what you had, and are easier for humans to read.

Records are simplified modules, and they're indexed exactly the same way: r:-a, r:-b. You're missing the second character of the indexing operator.

Unfortunately, convert(v, float) doesn't do anything when v is a string. To convert a string that represents a number to a float, you can use evalf(parse(v)). The whole code is

v:= GetProperty(TextAreaRadian,value);
exv:= evalf(parse(v));
SetProperty(TextAreaDegree, 'value', `if`(exv::numeric, sprintf("%6.2f", evalf(exv*180/Pi)), NaN))

This is my first time coding and embedded table, so there's no doubt a better way to do this.

The parse will give an error for strings that don't represent valid Maple expressions. Perhaps that's okay. If that's not okay, the error can be trapped with try.

This is called floating-point contagion. Once even one decimal point is placed into an expression, they tend to infect all expressions derived from that expression, and they tend to spread further around within those derived expressions.

Some tips:

1. Never use decimal numbers where a simple fraction will do. Most especially, don't do it in exponents.
2. Put off as long as possible the introduction of any decimal numbers at all into the computations. It's fine (and even a good idea) to enter required decimal data at the beginning of a worksheet and even to give those numbers names. Just don't use them in the early computations.

Explanation: Every time x passes through a layer of brackets, it counts as one function evaluation which removes one layer of quotes.

Suppose that I replace the square brackets with functions B1 and B2 (which'll be defined slightly differently from each other) and I use evalapply to define what is done with the (x) outside the functions. Try this code: 

restart:
B1:= (x::uneval)-> 'procname'(x):
B2:= x-> 'procname'(x):
`print/B1`:= `[]`:
`print/B2`:= `[]`:
`evalapply/B1`:= (B::uneval,  x::uneval)-> map(`?()`, B, x):
`evalapply/B2`:= (B,x)-> map(`?()`, B, x):
trace(B1, B2, `evalapply/B1`, `evalapply/B2`):
B1(B1(''z1+z2''))(x); lprint(%);
B2(B2(''z1+z2''))(x); lprint(%);

Does that help your understanding?

Status: It's a feature. I don't mean in the sense that it was explictly designed with this behavior in mind. I mean that this behavior is a natural and logical result of the overall design.

Both of these independence tests are quite easy computationally. The results under method= Pearson (the default) should be identical to those provided by Statistics:-ChiSquareIndependenceTest. I compute it explictly here for readers who just want to see the p-value computation without a lot of elaborate code for input error checking and output formatting. 

ChiSqIndTest:= proc(O::Matrix, {method::identical(Pearson, G):= 'Pearson'})
description "Returns p-value for Pearson's or G chi-squared independence test.";
option
   author= "Carl Love <carl.j.love@gmail.com>, 25-Oct-2018",
   reference= (
      "https://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test",
      "https://en.wikipedia.org/wiki/G-test"
   )
;
uses AT= ArrayTools, St= Statistics;
local 
   C:= AT:-AddAlongDimension(O,1), R:= AT:-AddAlongDimension(O,2), #column & row sums
   r:= numelems(R), c:= numelems(C), #count of rows & columns
   #matrix of expected values under null hypothesis (independence):
   T:= Matrix((r,c), (i,j)-> R[i]*C[j], datatype= float) / add(R)
;
   #p-value by either method:
   1 - St:-CDF(
      ChiSquare((r-1)*(c-1)), 
      add(`if`(method='G', 2*O*~ln~(O/~T), (O-T)^~2 /~ T))
   )
end proc:

Example of use:

DrugTrial:= <20, 11, 19; 4, 4, 17>:
ChiSqIndTest(DrugTrial, method= Pearson);
                       0.0402558477482379
ChiSqIndTest(DrugTrial, method= G);
                       0.0358413972184247

To get a numeric solution, you'll need numeric values for the parameters, an initial value for x (other than 0), and initial conditions. Then do

ode:= 
   diff(Y(x),x)+(A/x+B/x^2+C)*Y(x)+D4*Y(x)^2/x+E*Y(x)^3+F/x+G/x^2+H/x^3,
   diff(y(x),x) = Y(x)
:
ICs:= Y(1)=1, y(1)=1:
Sol:= dsolve({ode, ICs}, numeric, parameters= [A,B,C,D4,E,F,G,H]):
Sol(parameters= [1$8]):
plots:-odeplot(Sol, [x,y(x)], x= 1..2);

It's much easier than that. Simply do not load the Units package or any subpackage. Top-level Maple also understands units and the Unit(...) multiplier. By using Standard, you force the conversion to watts, because "hour" and "day" are not "standard" units.

I suppose that you simply want to treat the function as if it were a variable, without incorporating all the minutiae of the Physics package. Then it can be done in a one-line procedure using freeze and thaw like this:

Fundiff:= (e::algebraic, f::function)-> (F-> thaw(diff(subs(f= F, e), F)))(freeze(f)):

So, if you want the derivative of sin(f(t)) wrt f(t) to be cos(f(t)), then do

Fundiff(sin(f(t)), f(t)).

Applying this to the situation from your other thread, you could do

Fundiff~(SimplifiedMassMatrixDot, theta[2](t))

(The other thread is convoluted, with many questions having already been answered, so I agree with starting a new thread.)

I don't know if what I'm about to say is the cause of your problem, but even if it isn't, it seems likely to me that it will be the cause of a future problem.

I notice that you have made great use of the concatenation operator || to construct names (symbols in Maple-lingo). It is dangerous to use programmatically constructed symbols (constructed by any means, || or otherwise) explicitly in procedure code because All programatically constructed symbols are global. That includes symbols constructed by ||, cat, nprintf, parse, and convert(..., name). It doesn't matter if the "stem" of the symbol is declared local.

[Edit: I edited the above to highlight the distinction between names and symbols as those terms are used by Maple. Maple's symbols are defined at ?symbol (pretty much any sequence of characters, with enclosure in special quotes needed in some cases). Maple's names are a superset of symbols: names = symbols union indexed symbols.]

Here are two procedures for this. The first, SolveFunctionalRecurrence, computes the first N iterations of a functional recurrence, with N being an integer input. The second, PlotFunctionalRecurrence, simply plots all functions produced by the first.

SolveFunctionalRecurrence:= proc(
   #recurrence, solved for the highest-index function, with that highest index being just
   #a name:
   Rec::And(typefunc(name, indexed(name)), 1 &under nops) = algebraic,

   #parameter values (in form p::name= v::numeric) and initial functional conditions
   #in form (P::name[k::integer](t::name)= ...): 
   ICs::set({name, And(typefunc(name, indexed(integer)), 1 &under nops)}= algebraic),

   N::nonnegint #number of iterations to perform
)
local 
   n:= op(op(0, lhs(Rec))), j, k, F:= table(),
   ord:= -min(
      map(
         f-> op(op(0,f)), 
         indets(eval(rhs(Rec), n= 0), 'typefunc'(op([0,0], lhs(Rec))['integer'])) 
      )
   )
;
   for k to N do 
      F[k]:= eval['recurse'](Rec, ICs union {n= k, seq(F[k-j], j= 1..min(k-1,ord))})
   od;
   [entries(F, 'nolist', 'indexorder')]
end proc:

PlotFunctionalRecurrence:= proc(
   #output of SolveFunctionalRecurrence:
   F::list(And(typefunc(name, indexed(integer)), 1 &under nops) = algebraic),
   trange::name= range(realcons) #name and plot range of continuous variable
)
local k, n:= nops(F);
   plot(
      rhs~(F), trange, 
      'color'= [seq('COLOR'('HUE', (.85)*(k-1)/n), k= 1..n)],
      _rest
   )
end proc:   

And here's an example of it's usage on your functional recurrence. However, the results produced by your example initial function exp(-t) were too boring to be worth plotting, so I changed the initial function to the standard Normal distribution PDF, exp(-t^2/2)/sqrt(2*Pi).

F:= SolveFunctionalRecurrence(
   P[n](t) = (diff(P[n-1](t),t) + (lambda+mu)*P[n-1](t) - lambda*P[n-2](t))/mu,
   {P[-1](t)= 0, P[0](t)= exp(-t^2/2)/sqrt(2*Pi), lambda= 1, mu= 2},
   9 #number of iterations
);

PlotFunctionalRecurrence(F, t= -6..3, gridlines= false);

Note that my procedure requires that the recurrence be input in a form solved for P[n](t) with all functional terms on the right side having a lower index, and the number of initial functional conditions must match the order of the recurrence (even if that means that some must be specified as the 0 function, as is the case here).

The following one-line procedure removes n randomly selected elements from list L:

RemoveRand:= (L::list, n::nonnegint)-> (N-> L[[combinat:-randcomb(N, N-n)[]]])(nops(L)); 

There may be other ways that'll be overall more efficient if, say, L has thousands of elements or the operation is repeated hundreds of times.

Any plot can be animated, but you need to specify a parameter that varies with time. You have three potential parameters (r, theta, U), but you haven't quite specified which is supposed to vary with time. Does this code give what you want? I made U the parameter that varied with time.

plots:-animate(
   plots:-implicitplot,
   [
      U = r*cos(theta), r= 0..2, theta= 0..2*Pi, 
      coords= polar, axiscoordinates= polar
   ],
   U= -2..2
);