The keys to doing this are using the Units package, the function-definition operator ->, and the commands solve and evalf. A benefit of using Units is that if the units come out correct, that's a fairly good check that you've set up the problem correctly.

The additional key to doing it neatly is to use eval to substitute numeric values so that that step can be put off until the end of each computation.

Download airplane.mw

I've made an extensive investigation of this twice over the years. I haven't reinvestigated with Maple 2016, so this may have changed. What I found was:

• The primary order of layering of the objects when a PLOT structure is displayed is determined by the object type (CURVES, POLYGONS, POINTS, etc.).
• The order cannot be changed by changing the order that the objects appear in the PLOT structure.
• The objects most in the background are those of two topological dimensions such as area fillls.
• Next are objects of one topological dimension such as CURVES and the boundaries of POLYGONS.
• Next are objects of zero topological dimension: POINTS.
• Next is plotted text.
• The object most in the foreground are the axes.
• Within each class of objects, the order of layering can be changed by changing the order in the PLOT structure, which can in turn be controlled by the order that objects appear in a display command.

The plot that you want can be obtained through numeric differentiation with command fdiff:

fy:= y0-> fdiff((r-> rhs(r[-1]))@pds:-value(f(x,y), x=0), [1], [y0]):
plot(fy, 0..10, labels= [y, ``]);

You need to define an "epsilon cloud" around each integer so that if x is sufficiently close to an integer, then f(x) will be 1.

f:= x-> piecewise(abs((x-round(x))/x) < Float(6, 1-Digits), 1, 0);

The value Float(6, 1-Digits) = 6*10^(1-Digits) is often used as the target accuracy for floating-point computation.

Try any of

numtheory:-cfrac([1,1,2]);
seq(numtheory:-cfrac([1$k, 2]), k= 1..3);
Value(NumberTheory:-ContinuedFraction([1,1,2]));
seq(Value(NumberTheory:-ContinuedFraction([1$k, 2])), k= 1..3);

When using Maple's 2D Input, you can't have a space between a function name and the following left parenthesis. These extra spaces get translated as multiplication operators, which is why you see * after sphereplot and fieldplot above.

The trouble with your method is that you repeatedly add elements to existing lists or sets. Unfortunately, that operation is very inefficient in Maple. It is so natural to approach this problem recursively; however, I can't think of any way to efficiently do it recursively in Maple, even after thinking about it for hours.

Here is a way to do it using seq. This method is faster than Iterator. (In fairness, I must say that I learned this technique from a post by Joe Riel where he applied it to Cartesian products.) I also show a better way to do your indexing procedure TkSubSetList. And I also compare with VV's recursive method.

Download AllCombinations.mw

Indices, which are formed with square brackets, and different from subscripts, which are formed with double underscore. The two forms display pretty much the same in the GUI. However, a subscript is an integral part of the name to which it is attached, and hence it can't be substituted for an integer or anything else. A subscript is intended for display purposes, not programmatic purposes.

Additionally, you'll achieve better results with mul rather than product (lowercase p). Here's your worksheet, modified:

Download WhatHappensHere.mw

Here's a much-simplified version of your code:

funcs := [x, 3*x, x^2, 5*x];
for i to nops(funcs) do listM[i] := diff(funcs, [x $ (i-1)]) end do;
M := convert(convert(listM, list), matrix);
linalg:-det(M);

P.S. Here's a more-direct way:

M:= matrix([seq(diff(funcs, [x$(i-1)]), i= 1..nops(funcs))]);

You think that D[3](y)(L) refers to the 3rd derivative with respect to the first variable. That's incorrect: Rather, it refers to the 1st derivative with respect to the 3rd variable (regardless of the fact that y is a univariate function). There are several ways to represent the 3rd derivative with respect to the first (and only) variable:

(D@@3)(y)(L)
D[1,1,1](y)(L)
D[1$3](y)(L)

Likewise for the 2nd or any other order of derivative.

You need to spell out Heaviside and use exp(...) instead of e^.... The command is

inttrans:-laplace(Heaviside(t-1)*exp(3*t-3), t, s);

The Wronskian command is predefined. See ?Wronskian.

The linalg package has long been deprecated. It's functionality has been replaced by the packages VectorCalculus and LinearAlgebra.

The internal reference to h must be stored as z5:-h because h is in z5, so it needs to know where to look up h on re-evaluation of the expression. But you don't need to see the z5:- if you use a `print/h` procedure. Change z5 to this:

z5:= module()
option package;
export
     h:= F-> expand(evalindets(D(F)/F, specfunc(D), d-> op(d)*'h'(op(d))))
;
local
     ModuleLoad:= proc() :-`print/h`:= ()-> ':-h'(args) end proc,
     ModuleUnload:= ()-> unassign(':-`print/h`')
;
     ModuleLoad();
end module:

Of course, you can include any other locals or exports that you want in the module.

The key thing is that you can get any function f to prettyprint differently by defining a global procedure `print/f`. Since the procedure needs to be global, it needs to be assigned in a ModuleLoad procedure.

Symbolic antiderivatives in Maple are allowed to be discontinuous. This especially happens with complex functions. You can see this discontinuity in this plot (the red line is the real part):

F:= I*arctan(sqrt(exp(2*I*t)-1))-I*sqrt(exp(2*I*t)-1):
plot(([Re,Im])~(F), t= 0..Pi, discont);

The int command correctly adjusts for the discontinuity when it's a definite integral.

In the line that defines vNN, you have

vNN= ....

That needs to be

vNN:= ....

In the next line, use eval instead of algsubs:

pressure:= eval(nMom, vNN);