I just want to make a small terminology point, which has no affect on how this Question is Answered or coded in Maple: What you have here is not an "algorithm"; it's simply a recursively defined real-valued function. If you want to consider some method by which this function and others similar to it could be converted to a non-recursive form, or a method to reduce the total count of arithmetic operations involved in its evaluation (recursive or otherwise), then those would be algorithms.

@tomleslie Thanks for checking!

So, the string conversion can be done by

ExMat2:= (x-> `if`(x::numeric, x, sprintf("%A", x)))~(ExMat);

Note that the format is %A rather than the more-usual %a.

@mmcdara Vote up, also. I think that this is the most-important thing that you said (not to say that any other part was unimportant!):

• The problem comes from the fact that setting Y :=1/2*X1 + 1/2*X2  does not give Y the attributes of a RandomVariable.

See my Reply to VV. This crucial problem which you've pointed out (and is the crux of numerous other short-comings of Statistics) could be easily solved in modern Maple if RandomVariables were implemented as objects and the arithmetic operators (such as * and + in the above expression) were overloaded with respect to that object. At the moment, I think that the majority of the work to code these objects could be done without touching the underlying Statistics package, which makes this work relatively easy.

@vv Vote up. However, I don't totally agree with your statement "Statistics is a modern package". Maple 16 (some 9(?) years ago) refined the concept of modules to include objects. RandomVariables should be implemented as objects, and Statistics (which is a heavily used package) desperately needs to be updated along these lines. 

@subzero I'm sorry to learn that it doesn't work for you. What Maple version are you using?  Although I'm not aware of using any of the numerous relatively new Maple features in that code, someone else here might recognize something. 

Please run this diagnostic: After entering the procedure code, give the command

trace(specdifforder);

Then define the expression that has the derivatives, and run at least one case that works and at least one that doesn't. There'll be a lot of hard-to-understand output (which'll nonetheless be very useful for me). Save the worksheet and post it here as an attached file. Files can be uploaded and attached by using the green uparrow on the toolbar. If you get an error message during the upload, just ignore it and proceed; it's likely that the attachment will work anyway.

@mmcdara You asked, rhetorically:

• But the things are not that simple...
  As changing 0.5 to 1/2 "works", one can ask ourselves  if, for instance, a*X1 + (1-a)*X2 still "works" when a is a rational between 0 and 1 (infering this doesn't work if a is a floating point number between 0 and 1)?

In the case of a any exactly represented number in 0..1 (whether rational or irrational), it is indeed "that simple": Maple will correctly compute PDF(a*X1 + (1-a)*X2, t) and the floating-point evaluators used by plot will be able to achieve the required precision over the whole interval.

@vs140580 So what makes this output any different from an ordinary Maple matrix? Why won't you give an example worksheet?

Okay, from your other Questions, I'm guessing that Python output always appears in Maple as a string? Is that right? Or at least it's always a string the way that you do it. So, do you simply want to write one of these strings, which represents a matrix, to a text file? Or do you want to convert it to a Maple matrix and then write that to a file? 

@P2prod You wrote

• So the above expression is not a function but a value depending on x and y (piecewise). Isn't it?

No, quite the opposite. The expression does not depend on x or y, and the fact that piecewise is involved is irrelevant.

• If I write evalf(A) I get 2.956602889

Try that again, starting from a restart. The exact value is 0, and evalf(A) returns 0. for me. Plotting the integrand may make that obvious:

plot3d(
    2*piecewise(x < 1, x, 1 < x, 2 - x) * 
        piecewise(y < 3/2, y, 3/2 < y, 3 - y) *
        sin(Pi*x/2)*sin((2*Pi*y)/3), 
    x= 0..2, y= 0..3
);

Do you also find the following surprising?

int(piecewise(x<0, x, x > 0, -x), x= -1..1);
                               -1
int(x, x= -1..1);
                               0
int(-x, x= -1..1);
                               0

It seems to me to be totally analogous to your double integral.

Other than the trouble with evalf(A) (hopefully you'll get 0 next try), the issues that you're struggling to understand are 100% mathematical. I will be happy to explain those mathematical issues. Start by plotting the following 3 functions on the interval -1..1:

f1:= x:
f2:= -x:
f_pw:= piecewise(x < 0, x, x > 0, -x):

Looking at each plot, determine (purely from geometric considerations) what its integral from -1 to 1 should be. I have a feeling that just doing this exercise will cause "a lightbulb to go off in your head", and you'll understand the whole issue. Either way, let me know.

Please provide a Maple worksheet that produces and contains an example of the matrix. In light of your other Questions today, I don't yet see what. if anything at all, makes these any different from ordinary Maple matrices.

@P2prod Using variable upper limits of integration, I get this:

plot3d(
    (2*Int(3*piecewise(x < 1, x, 1 < x, 2 - x) * 
           piecewise(y < 3/2, y, 3/2 < y, 3 - y) *
           sin(Pi*x/2)*sin((2*Pi*y)/3), 
           [y = 0 .. a, x = 0 .. b])
     )/3, b = 0 .. 2, a = 0 .. 3
);

@roethlisalayna 

If A is the set of odd numbers, then

select(isprime, A)

@vs140580 You didn't answer my second question: Does trunc~(A) do what you want?? If that doesn't help, how about

trunc~(rtable(A, datatype= anything));

@AmirHosein Sadeghimanesh You wrote:

• Recently a colleague showed me that it is better to use `__i` instead of `[i]` for indexed variable names.

I wouldn't categorically say that one is better than the other: They each have their uses. One thing that you can never do is make independent assignments to both x and x[i], but you can do it with x and x__i.

Your cat(j, '__', i) can be replaced by cat(j__, i). One drawback to using cat is that the variables so constructed are always global. This is true of all programmatically constructed names, regardless of whether the construction is done with cat, nprintf, ||, or parse.

@lcz I certainly expected some comments on this procedure.