Since you are a student of number theory, you should immediately prove the following three elementary propositions about divisors, which'll lead to one efficient method of solving your posted problem (efficient for the vast, vast majority of practical cases). I both came up with and proved these propositions as my first formal written proof as a 16-year-old college freshman.

1. The number of divisors of a positive integer that are less than its square root exactly equals the number that are greater.

2. The number of divisors of a positive integer n is odd iff n is a perfect square.

3. If a positive integer n has no divisors d such that 1 < d <= sqrt(n), then n is prime.

(The only criticism that I received on that freshman work was due to the fact that in my naivete, I called them "factors" rather than "divisors".)

I'll let that sit for a few hours, then I'll post two algorithms, one based on sort that uses proposition 1 above and one that uses min and select, akin to Mr Heinz's but much faster. Now, based on my theoretic knowledge of algorithmic asymptotic time complexity, I know that there's guaranteed to be some unimaginably HUGE number n such that the min/select method is faster than the sort method, but I'd guess that the smallest such n has millions of divisors or more. That why I said "practical cases" above.

I've read the paper that you linked the PDF of. The ODE system in the paper is a two-point BVP of four second-order ODEs. Maple's dsolve comes with four solvers dedicated to two-points BVPs. RK4 is an IVP solver (and, as I said above, it's been long ago superseded by several better IVP solvers, the most notable being RKF45). To use an IVP solver to solve a BVP, one must construct an elaborate "shooting method":

1. let some of the initial conditions be chosen arbitrarily
2. solve an IVP
3. evaluate the solution at the other boundary point. If they're close enough, you're done; if not, refine the initial conditions and go back to step 2.

The authors only give the barest hint of this process by mentioning Newton's method. Since BVP solvers are much flakier than IVP solvers, the above process is occasionally necessary. In the case of this BVP, it most definitely is not. So, it seems to me that in lieu of any explanation for why they used RK4 + shooting + Newton---and they give none---the authors are foolish to use this elaborate method. They should provide you with their Maple code.

In the following worksheet, I easily and systematically generate the first 56 plots from the paper, that being all the plots in figures 2 - 5. I could've made the rest of the plots, but I got tired of it, and I think that you can and should follow my example, which'll easily get you the remaining umpteen 2D plots. If you care to generate the 3D plots also, I can help you with that just as systematically.

Note that my plots are much smoother than those shown in the paper, especially when you view them in Maple rather than on MaplePrimes. That's because the BVP solver is automatically choosing the interval width needed to make them smooth. It's all automatic! All the default dsolve settings worked perfectly (although that's often not the case with BVPs).

This BVP is homogenous (i.e., the independent variable does not appear) and nearly linear. The only nonlinearities are two quadratic terms neither of which involve the highest-order derivatives. That's about as easy as it gets for practical numeric BVPs that we see on MaplePrimes.
 

Download nanofluid_BVP.mw

[Before you download the above, note that there's an updated version in a Reply below. The update includes everything that's in the original.]

Vote up for a great Question.

The symbols &under and `&+` are defined on the page ?type,structured, which is one of deepest pages in whole help system and one of my most commonly reread (the other being ?operators,precedence). I must add that both &-operators and their use in structured types has been around for a long time---at least 18 years, which covers the extent of my usage.

The & is being used in two different ways above. Let me first explain &under. Any symbol (with a few exceptions to disallow ambiguity or excessive weirdness) that begins with & can be used as a binary infix operator, regardless of whether the symbol's meaning has been predefined. If I type A &foo B, then the function `&foo`(A,B) is constructed. If a procedure `&foo` has been defined, then it is called; otherwise the construction remains an unevaluated function like any other. Note that the quotes are required for the prefix form `&foo`(A,B) and that they are not used in the infix form A &foo B, and these two expressions mean exactly the same thing.

So, the call type(e, T &under F) becomes type(e, `&under`(T,F))---a parameterized type. Like any old-style (i.e., pre-TypeTools) parameterized type that doesn't have a built-in handler, this invokes `type/&under`(e, T, F). Use showstat to see the one-line code of `type/&under`. Then you'll see that the original type command is equivalent to type(F(e), T). So why would one use the former? Because now T &under F is itself a type that can be used with the numerous fundamental commands that require a type argument, such as indets, hastype, membertype, procedure parameter restrictions, etc.

(Typing up explanation for `&+` now. Need to post so the above doesn't get lost. One wrong keypress in MaplePrimes and your whole post is gone...

Okay, here it is....)

Note very carefully that & is not itself an operator! It's just a character and has no semantics by itself. Any characters---no matter how weird (including control characters, spaces, escape characters, unicode, quotation marks, etc.)---can be part of a Maple symbol as long as the symbol is enclosed by ` `. I've worked on a large system where all the internally generated symbols (mostly dummy variables of integration) were Chinese characters. If the symbol's first character is &, and the rest are normal printing characters, then the symbol can used infix (there are also some ambiguous constructions that are not allowed).

The purpose of `&+` is to have something that looks kind of like `+` but has no procedure assigned to it. If you were to instead use `+`(...), then the addition would immediately evaluate before the type command saw it. The command type(e, `&+`(T1, T2)) returns true iff e is the addition of something whose first operand is type T1 and whose second operand is type T2 (and they must be in that order). The type `&+` is built-in to the type command; you can't see its code.

The type system is the essence of symbolic programming.

I think that "unknown" is a poor and confusing choice of words for this. I think that the term anonymous procedure is the most-used term for this concept in computer science, and even in the Maple documentation.

Of course the system knows what procedure produced the error. The issue is that in Maple it's easy and common to create and use a procedure without ever giving it a name. It's only the procedure's name that's unknown, not the procedure itself.

For example:

map(x-> sin(1/x), [0, 1]);

Here, the procedure x-> sin(1/x) is anonymous, and the error message is like what you show.

To make matters worse, Maple actually gives names to these anonymous procedures. They're all of the form `unknown/H` where H is some string of hexadecimal digits. If this name were shown in the error message, it would make things so much easier! Then you'd be able to view the errant procedure by applying eval to its name. I'm try to figure out a way to get the name. I know it'll involve using debugopts(callstack), but I haven't worked out the details yet.

@Gabriel samaila

The first argument to your dsolve command should be sys1, not dsys. Otherwise, there was no point in creating sys1, right? As the warning message says, you can also use the parameters option. That would let you specify the value of the parameters after running the dsolve command.

Why do you want to use Runge-Kutta 4 specifically? There are newer, better methods. These newer methods select and vary the stepsize automatically, basing their choices on getting the absolute and relative errors below certain user-specifed tolerances.

Any number of plots of the same dimensionality can be combined with plots:-display, for example:

plots:-display(
   [
      Statistics:-Histogram(
         Statistics:-Sample(Normal(0,1), 9^5),
         color= aquamarine
      ),
      plot(exp(-z^2/2)/sqrt(2*Pi), z= -4..4, thickness= 5)
   ],
   title= "       The standard normal curve"
);

Notice that options that apply to individual plots (such as thickness= 5) go with the individual plot commands, and options that you want to apply globally (such as title) come after the individual plot commands.

Any Maple plotting command creates a data structure that can be assigned to a variable or manipulated otherwise. So the above could also be done like this:

P1:= Statistics:-Histogram(Statistics:-Sample(Normal(0,1), 9^5), color= aquamarine):
P2:= plot(exp(-z^2/2)/sqrt(2*Pi), z= -4..4, thickness= 5):
plots:-display([P1,P2], title= "       The standard normal curve");

These things are better done in Maple without using for. Anyway, we've shown you several for loops by now, so I think that you should be able to come up with something on your own. Here's my solution for the first problem:

Insert:= (x, pos::And(posint, satisfies(n-> n < nops(L)+2)), L::list)->
   [L[..pos-1][], x, L[pos..][]]
:

If you're studying computer science, you should learn that changing lists is inherently inefficient when you need to shift large numbers of list elements to new positions. And if you're studying Maple specifically, you should learn that it's especially inefficient in Maple, and that Maple offers other indexable order-preserving container structures (such as Vectors and tables) that are more efficient when those structures need frequent updates.

The answer is that you simply use F inside P. I removed the evalf, which only complicates things and makes them more difficult to debug. Other than the evalf, your code is equivalent (unless I messed something up) to this:

lambda:= .2:                              
mu:= .3:
                              
F:= n-> exp(-lambda*t)*(lambda*t)^n/n!:

P:= proc(n) 
local ret:= F(0), i; 
   for i to n do
      if i = 1 then ret:= ret + (diff(F(0),t)+lambda*F(0))/mu
      elif i < n then ret:= ret + (diff(F(n-1), t) + (lambda+mu)*F(n-1) - lambda*F(n-2))/mu
      else ret:= ret - (diff(F(n),t) - lambda*F(n-1))/mu
      end if
   end do;
   ret
end proc:

To test that, it'll help if you don't assign the values of lambda and mu. Let them remain symbolic.

But, I think that you made a mistake in P. It bothers me that the code after elif i < n does not depend on i. I think that you meant in that line F(i-1) and F(i-2). I can incorporate that change into a further simplified procedure:

restart:
#Leave lambda and mu unassigned for testing.
#lambda:= .2:                              
#mu:= .3: 
                             
F:= n-> exp(-lambda*t)*(lambda*t)^n/n!:

P2:= (n::nonnegint)-> 
  (`if`(n > 0, diff(F(0),t)+lambda*F(0), 0) +
   add(diff(F(i),t) + (lambda+mu)*F(i) - lambda*F(i-1), i= 1..n-2) -
   `if`(n > 1, diff(F(n),t) - lambda*F(n-1), 0)
   )/mu +
   F(0)
: 

Please test that. P2 is my replacement for your P.

My favorite way to do this uses irem. If n and d are positive integers, then irem(n, d, 'q') returns the remainder from dividing n by d and stores the quotient in q. So, my procedure is

Bin:= proc(n::posint)
local q:= n, bit, k, i;
   for k while q <> 0 do
      bit[k]:= irem(q, 2, 'q')
   od;
   [seq(bit[i], i= 1..k-1)]
end proc: 

The order that the bits are returned by this may seem backwards to you; however, this is the order that's most convenient for most subsequent computational uses. If you want the other order, simply change bit[i] to bit[k-i] in the last line of the procedure.

Here's a simple way to construct your matrices, multiply them, and force through the explicit multiplication of parenthesized factors:

n:= 6: # order of matrix
# exp(I*z) = cos(z) + I*sin(z)
A:= Matrix(n$2, (i,j)-> exp((1-i)*(j-1)*2*Pi*I/n)); 
expand~(A.A^*);

Notes (some of which were covered by others above):

• # indicates the beginning of a comment, which continues to the end of the line.
• sqrt(-1) in Maple is I, not i, unless you change that default with interface(imaginaryunit= ...).
• For a matrix A, the conjugate transpose is A^*. Since your matrix is symmetric, it's effectively just the conjugate.
• The expand command forces through the multiplication of parenthesized factors.
• For any operation OP that you want to apply to the elements of a matrix M, you can make the operation apply to all the elements by OP~(M).
• One of the most important identities in all of mathematics is exp(I*z) = cos(z) + I*sin(z). This doesn't change your matrix at all; it merely simplifies how you input it.

You can upload a worksheet to this website by using the green up-arrow on the tool bar of the editor. Nearly all of us who answer questions here prefer this to (or in addition to) screenshots because that way we don't need to retype your input into Maple.

In the general case that the Vs are arbitrary functions (expressions in a single-variable x, or constants), you can use piecewise as in

plot(piecewise(seq([x < L[k], V[k]][], k= 1..numelems(L))), x= 0..L[-1]);

If the Vs are just line segments or just constants (as in your example), then a more-efficient solution is possible, if you care, although the piecewise solution is not notably inefficient.

Note that I'm simply using the plot command; plots[multiple] isn't needed, or desired, for this.

Of course, my command above assumes that

1. L and V are indexable structures with the same number of elements and indices starting at 1 (such as Vectors or lists);
2. the Ls are distinct positive constants sorted in increasing order;
3. you wish the domain to start at x=0.

If any of those are not true, my command can be adjusted to accomodate that.

In addition to what Acer said, which is totally correct, I thought that I might give a more-thorough answer to the "How?" question that you directly posed at the end of your Question. Take a look at this simple recursive code, which is what's ultimately executed by your convert command:

showstat(`convert/list`::convert_rtable_to_nested_list);

You'll see that there's no magic efficiencies, no special provisions for sparsity, no special provisions for order, and no special handling of zeros. It's just a recursive loop whose base case is simple indexing, as in A[...].

Integration with respect to functions as opposed to with respect to simple variables is not something handled by int (well, at least not in its common usage). If you assign expressions to x and y, then you can't integrate with respect to x or y.

Most packages are implemented as modules. Anything that can be accessed in the form A:-B comes from a module named A. Modules have two types of local variables: regular (private) locals and (public) exports. Ordinarily, only exports can be accessed with the :- operator. But, if you issue

kernelopts(opaquemodules= false);

then you'll be able to access private locals also.

Another way to see module A's local procedure B is

showstat(A::B);

This was only implemented a few versions ago, so I don't know if you have it. Regardless, the opaquemodules trick will definitely work.

Note that it's ithratB, not ithratb. Also, showstat(ithratB) will never work. It either has to be showstat(numtheory:-ithratB) or showstat(numtheory::ithratB) (or <*cringe*> showstat(numtheory[ithratB]) (*footnote)). Dropping the module-name prefix is only allowed for exports and then only in the context of a with or use command or a uses clause.

(*footnote): The use of the syntax A[B] as a general substitute for A:-B is something that I strongly disdain. The A:-B is precise: It can only mean that B is literally and explicitly, without any further evaluation, a local or export of module A. The syntax A[B] can mean a great many things and even if it's known that A is a module, it's B's evaluated value that gets looked up. The A[B] is only acceptable to me when B needs to be evaluated first. Unfortunately the Maple documentation is loaded with the inappropriate use of this syntax.

Something that I occasionally find useful along these lines is kernelopts(memusage). It returns a listlist (which, oddly, displays as a Matrix---I guess kernelopts has special privileges of display?) of three columns: The first is labels for 63 categories of possible objects stored in memory, the second is the number of objects stored in each of those categories, and the third is the number of bytes stored in each of those categories. So, you can call it before and after calling a suspect procedure and take the difference of the two readings.

I believe---but I'm not sure---that this command gives you a true partition of all stored data, i.e., every byte of data that you currently have stored appears in exactly one of the 63 categories.