@Spartan 

Not every module local or export is a procedure. Some of them are just unassigned names, just like the locals of most procedures. Or they can be any other Maple data structure. In the pristine state (i.e., just after a restart), currEqn is just an unassigned name. I suspect that if you recheck after actually running DiffTutor, it'll be assigned the expression that you were trying to differentiate.

I had intended that the procedure that I gave would actually print out all the procedures, not just list their names. But due to a dumb mistake, I wasn't in the pristine state when I tested. Here is what I intended for my procedure:

ModulePrint:= proc(M::`module`)
uses ST= StringTools;
local Ns, N;
     interface(verboseproc= 2);
     kernelopts(opaquemodules= false);
     try unwith(M) catch: end try;
     Ns:= eval~([op(3, eval(M))]);
     Ns:= ST:-Split(sprintf("%q", op(3, eval(M))), ",");
     for N in Ns do
          lprint(N);
          parse(sprintf("print(eval(%s))", N), 'statement')
     end do
end proc:

ModulePrint(Student:-Calculus1:-DiffTutor);

This'll list all the procedure code for you; you should need no other commands. If it doesn't, let me know.

@bpourhamzeh I doubt that I would've figured that out simply by looking at the menus.

If you have a restart command at the top of your worksheet (which is almost always a good idea), keep it in its own execution group.

@ecterrab 

Have you looked at the examples in my previous Reply? How do you explain them?

No, my question is not yet clear to you. I do not want the indices to start at different values! That would be ridiculous. But I am trying to show you that your system as currently implemented DOES cause the indices to start at different values even when that is not what the user intended. Did you think that the examples that I posted in my previous Reply represented how I wished the system to behave? No, they are examples of actual execution, of how it DOES behave.

What I want, if possible, is for both indices to start at 1. If they must start at 0, I can work with that.  But the behavior shown in my second example must be a bug, right?

You said:

where the lower limits are always determined by the value of the right-hand-side of the multi-index sum

This has me baffled. Did you mean the left-hand side? My first example above clearly shows that making the left-hand side 1 does NOT make the indices start at 1.

Haley,

I am very disturbed about not getting a reply from you. I put a lot of work into my Answer, and I think that it is generally useful. Does it answer your question?

@tomleslie You asked:

Maybe if someone could tell me for which type of argument my simple-minded overload gives a different answer from Carl's

Those types of arguments are named tables, named modules, and named procedures. This is because of the "last name evaluation" rule that applies to those (see ?last_name_eval). I could've also coded it as

proc(C::And(name, Not(last_name_eval)))...

I hope that this explains the &under concept once and for all: type(X, T &under P) if and only if type(P(X), T). In the case at hand, I don't simply want to require that C be of type name; I want to require that eval(C) be of type name.

Your "simple-minded" attempt was very good. It would've been thwarted by a named table argument, but this last name evaluation rule is a very special case.

You've been asking very similar Questions for a few weeks, as if you're not learning from the Answers that you've been given. So, what's going on? You should post your attempted Maple work. It helps us to see what you're doing wrong.

Excellent Question.

What is the density of the matrices (the proportion of nonzero entries)? Is there a highly structured pattern to the nonzero entries, such as tridiagonal? What is the average and maximal degree of the polynomials for the dimension = 674 case? Are the polynomials themselves sparse? If so, what is their average density? What is an upper bound for the integer coefficients? Does your algorithm require multiplying different such matrices or just raising matrices to powers?

Time the run of your 34-second example like this

CodeTools:-Usage(LinearAlgebra[MatrixMatrixMultiply](A,A)):

and report the results. (I'm particularly interested in the garbage collection time.) It doesn't need to be exactly the same example that you used earlier.

@al-faraj 

Your function shows the position of the particle, not the pace. To get the pace, you need the derivative of the function.

solve(diff(5*t^(1/2), t) = 2, t);

@tomleslie It won't work because the parser objects to having a loop (or any statement) split across multiple execution groups.

@acer 

This probably already occurred to you, but you can avoid the overhead of calls to an anonymous procedure v-> F.v by using map2. Change map( v -> F . v, CCXYZ_D50 ) to

map2(`.`, F, CCXYZ_D50 )

@Markiyan Hirnyk Thank you very much for the encouragement. That's a beautiful painting. The clouds are similar to the Gaussian Random Fields with large convolution radii. I find them also beautiful. Here's an example of a GRF(512, 128) (the 128 is the convolution radius). (Sorry for the excessive whitespace---these are 3D plots projected onto one of the coordinate planes; I don't know how to eliminate it.) It looks like this could be a model for how the shapes of continents, islands, and oceans develop.

In addition to the bracket problem, you need to change i. In Maple, sqrt(-1) is represented by default as I. You can change that to i if you want, but it requires an explicit command.

@Markiyan Hirnyk 

Okay, how about this?

restart:

Digits:= 15:

GRF:= proc(n,r)
uses AT= ArrayTools, IT= ImageTools;
local
     GRF:= IT:-Convolution(
          AT:-RandomArray(n, 'distribution'= 'normal'),
          Matrix(2*r+1$2, (i,j)-> `if`((i-r-1)^2 + (j-r-1)^2 <= r^2, 1, 0))        
     ),
     Pos:= map[evalhf](`>=`, GRF, 0),
     Neg:= map[evalhf](`<`, GRF, 0)
;
     (IT:-Preview@IT:-HSVtoRGB)(
          IT:-CombineLayers(.6*360*Pos, IT:-FitIntensity(Pos*~GRF), Pos) +~    
          IT:-CombineLayers(.1*360*Neg, IT:-FitIntensity(-Neg*~GRF), Neg)
     )
end proc:

GRF(200, 13);

The underlying greyscale image is the same as what I posted previously; only the color encoding is different. The two colors can be easily changed to any pair on the Hue (the H of HSV) scale, i.e., colors that correspond to specific wavelengths of light. I could put a gamma encoding on the Saturation (the S of HSV) layer to raise some of the whitish tints. (The normally distributed values close to zero are the colors close to white.) Make sure to view the image from a variety of oblique angles. It can make a big difference because of how the RGB layers are laid out in your monitor screen.

Here's a more-sophisticated color encoding that is fairly close to the Mathematica-generated plot that you posted:

GRF:= proc(n,r)
uses AT= ArrayTools, IT= ImageTools;
local
     GRF:= IT:-Convolution(
          AT:-RandomArray(n, 'distribution'= 'normal'),
          Matrix(2*r+1$2, (i,j)-> `if`((i-r-1)^2 + (j-r-1)^2 <= r^2, 1, 0))        
     ),
     Pos:= map[evalhf](`>=`, GRF, 0),
     Neg:= map[evalhf](`<`, GRF, 0)
;
     (IT:-Preview@IT:-HSVtoRGB)(
          IT:-CombineLayers(
               .6*360*Pos*~IT:-FitIntensity(Pos*~GRF, .85..1.2),
               IT:-FitIntensity(IT:-FitIntensity(Pos*~GRF)^~.75, .08..1),
               .93*Pos
          ) +~    
          IT:-CombineLayers(
               .1*360*Neg*~IT:-FitIntensity(Neg*~GRF, .4..1.5),
               IT:-FitIntensity(IT:-FitIntensity(-Neg*~GRF)^~.75, .08..1),
              .93*Neg
          )
     )
end proc:

GRF(200, 12);

This uses two ranges of Hues (positive (blue): 0.51..0.72, negative (orange): 0.04..0.15; on a 0..1 scale) rather than two specific Hues. A gamma (or exponent) of .75 is applied to the Saturation layer and then the minimal value is set to .08 (on a 0..1 scale). The Value layer is constant .93 (on a 0..1 scale).

Again, this is just a translation (plus color encoding) of the first of the many Matlab programs given in that paper, and this is not the algorithm for the plot that you posted, which, I believe, was done with FFT rather than convolution. I may translate some of the others.

@tomleslie I really like that as a quick and dirty solution, but it's difficult to generalize it to an arbitrary number of matrices.

@Carl Love 

So, how did I get those expressions for the partial derivatives of F, the key to the previous computation? Working with multivariate chain rules and the inverse function theorem to get derivatives by hand is difficult. Here are two ways to do the symbolics with Maple:

F:= (y,z)-> RootOf(MF(x,y,z)-11, x):
D[1](F)(y,z);

D[2](F)(y,z);

The following way is perhaps easier to understand:

restart:
alias(F(y,z)= RootOf(MF(x,y,z)=11, x)):
diff(F(y,z), y);

and likewise for diff(F(y,z), z).