You need to substitute a procedure (colloquially called a function) for x[1] rather than substituting an expression for x[1](whatever). Like this:

eval(junk, x[1]= ((a,b,c)-> R(b,c)*sin(omega*a + phi(b,c))));

Always use eval for these high-level calculus substitutions. Using subs isn't necessarily wrong, but eval is safer. The use of subs should be reserved for low-level purely syntactic substitutions.

I am highly skeptical of Kitonum's Answer. I don't doubt the computed critical value of D(theta)(0)---I got the same value. But I do doubt the solution curves. And I doubt the existence of a finite T such that theta(T) = Pi and D(theta)(T) = 0. My intuition is that diff(theta(t), t) (for t > 0) is a strictly positive, strictly decreasing function whose limit at infinity is 0.

Here are some computations to support these hypotheses:

Download Highly_unstable_ODE.mw

In general, there's no way to do it. But many builtin (or kernel) functions come from open-source libraries. For example, reading the code for isprime, one sees that it uses (for some inputs) gmp_isprime. That's part of the GMP open-source library. So, if you Google "gmp isprime", you'll quickly find that code in C.

How about saving everything to one .mla file?

save H||(1..n), "H.mla";

I just read your other Question, regarding Grid, and I now realize that you probably really need nine separate files. That can be done for the kth file by

save cat(H,k), sprintf("H%d.mla", k);

By the way, you may be confusing things by calling these .mla files. Use another extension. If you use .m, the files will be compressed into Maple's own text-based but not-human-readable format, which'll probably save an iota of time.

That needs to use nested seqs:

seq(seq([A[i],B[j]], i= 1..5), j= 1..4);

or

seq(seq([A[i],B[j]], j= 1..4), i= 1..5);

A better way, which doesn't require hardcoding the list sizes, is

seq(seq([i,j], j= B), i= A);

A method that works for an arbitrary list of lists L, each of arbitrary length, is

L:= [A,B]:
CP:= combinat:-cartprod(L):
'CP[nextvalue]()' $ mul(nops~(L));

You can't assign to a procedure's parameter, n in this case. You need to copy the parameter to a local variable, and then assign to that.  Like this:

cs:= proc(N::posint)::nonnegint;
description "Find the number of steps for a Collatz sequence to reach 1.";
local n::posint:= N, q::posint, r::identical(0,1), count::nonnegint;
     for count from 0 while n > 1 do
          q:= iquo(n, 2, 'r');
          n:= `if`(r=0, q, 3*n+1)
     end do;
     count
end proc:

Also, your count was off by one to my mind: cs(1) should be 0.

The properties (in the strict Maple sense of "properties") of an object can be inquired with the about command; however, this is only useful if you've used the assume command, which is the command which gives objects properties.

Somewhat similar to properties are attributes, which can be inquired with attributes. An object will only have attributes if you give it attributes with setattribute. You can give attributes whatever meaning you want.

To see all subexpressions of an expression e, use indets(e, anything). The anything can be replaced with any type expression (in the strict Maple sense of "type" (see ?type)). For example, all the integrals in e can be found with indets(e, specfunc({Int,int})). This makes indets one of the most useful and most commonly used commands.

To have e broken down into a logical (but not visual) tree which includes the smallest detail of information, use ToInert(e). This is often too much information to wade through by reading; however, this information can be manipulated programmatically. Some of this information may be arranged in a logical fashion that doesn't exactly correspond to the internal representation of the object, but does correspond to the logical representation of the object.

To have e broken down into a more-visual tree, but which uses some more-cryptic abbreviations, use dismantle(e). This information, as far as I know, corresponds precisely to the internal representation. This information is quite difficult to manipulate programmatically.

Ultimately, everything in a computer is stored as numbers. If you really want to go down to that level, use disassemble(addressof(e)). This will break it down one tree level only. Each given number can be broken down further with disassemble if it's an address. It's a good idea to follow the dismantle tree when you are doing this.

You should put your Maple version in the header of your Question. There's check-off boxes and a pull-down menu for it.

Your zero-one Matrix can be created like this:

n:= 4:
subsMatrix:= <seq(x[], x= Iterator:-CartesianProduct([0,1]$n))>;

Your list of substitution equations can be created directly from the iterator like this:

A:= [seq(a[i], i= 1..n)]:
subsEqs:= [seq(
     select(e-> rhs(e)=0, A=~ convert(x[], list)),
     x= Iterator:-CartesianProduct([0,1]$n)
)];

So note that there's no need to create the subsMatrix if its only purpose was to be an intermediary for the creation of the subsEqs.

Your error? It's probably that you need to say LinearAlgebra:-RowDimension rather than simply RowDimension and likewise for ColumnDimension.

The Iterator package comes with Maple 2016. If you have older Maple, let me know.

In bcs, you need to put parentheses around f1 in the first boundary condition and around f3 in the fourth boundary condition.

The command indets(sys2 union {bcs11} union bcs, name) can be used to find such problems. It this case, it told me that D was being considered a name (as well as an operator).

You asked:

How can I see the order of vertices of square? Because, perhaps you sorted....

Kitonum's code maintains the four vertices of each square as a set, and all the squares themselves are in a set. Sets in Maple are sorted into stored in an order determined by Maple (mostly lexicographic order). To maintain the order of the vertices within the squares, it's necessary to store the squares as lists instead of sets. This causes duplicate squares, so they need to be divided into equivalence classes by their equality as sets, and have a representative of each equivalence class selected. A command that can form equivalence classes is ListTools:-Classify.

Here's a procedure that does this for any sphere with integer center and radius. This is based on Kitonum's code but is generalised to work for any sphere with integer center coordinates and radius. It takes significantly more time because of the equivalence class finding. The coding is more complex due to the many subtle list-set-table conversions required.

Download AllDistinctIntegerSquares.mw

I can get Maple to put it into the following relatively simple Sum form, which is not as simple as the one shown in the paper. Note that the only difference in the input is that I used Sum with an uppercase S. That's the inert form of sum, but diff still knows what to do with it.

The following commands have inert forms that the other commands know how to handle. In all cases, the inert form is the same word capitalized: sum, product, diff, int, eval, limit.

Download diff_double_sum.mw

Note that the formulae that you entered aren't quite as general as those in the paper. In the paper, the sums have coefficients that depend on n and m.

Is this what you're looking for?

linspace:= proc(r::range(complexcons), n::posint)
local a:= op(1,r), h:= (op(2,r)-a)/n, k;
     [seq(a+k*h, k= 0..n)]
end proc:

linspace(0..1., 4);

     [0., .250000000000000, .500000000000000, .750000000000000, 1.00000000000000]

The name linspace comes from a very similar function in Matlab.

Why not simply have a variable, Nlists, that you add 1 to every time you find a list?

Nlists:= Nlists + 1;

Since your task seems so large (you mentioned 11 trillion lists in another Questions), you should describe it in more detail so that we can discuss its feasibility. Perhaps there is already a combinatorial formula to count it. Perhaps your lists could be generated by the Iterator or combstructs package. Note that you can't fit 11 trillion lists on even a large disk. You'll need to use an iterator of some sort. That is, the lists need to be created, examined, counted, and garbage collected in a loop. Maple will handle the garbage collection on its own.

The bottom of your Question is about the data structure that you should use. I think that a DataFrame would be good for this. A DataFrame is like a Matrix except that the rows and columns are indexed by labels rather than by numbers. I think that your columns should be labeled by weather phenomena and the rows should be labeled by date and hour. See ?DataFrame.

There are no functional differences. If there were, we'd occasionally have to ask people when they asked Questions here what edition they had so that we could answer. That has never happened as far as I'm aware. In Maple V (about 15 years ago), there was a small difference in the student edition: Precision was limited to 100 digits. I'm not aware of any differences after that.