You should define the Matrix with symbolic a (i.e., before assigning to a).

restart:
M:= <a, 2*a; 3*a, a^2>;

a:= 5:
M;

unassign('a');
M;

Is the plot below what you want? Note that I've added a third ODE to the system. This specifies a function H for the integral of h. I've also added a corresponding boundary condition, H(0)=0. This process is faster and more accurate than direct numeric integration of h.

eq1:= diff(f(eta),eta$3)+f(eta)*diff(f(eta),eta$2)-diff(f(eta),eta)^2+1:
eq2:= diff(h(eta),eta$2)+f(eta)*diff(h(eta),eta)-diff(f(eta),eta)*h(eta):
eq3:= diff(H(eta),eta) = h(eta):
bc:= f(0)=0, D(f)(0)=0, D(f)(6)=1, h(0)=0, D(h)(6)=1, H(0)=0:
A1:= dsolve(
     {eq||(1..3), bc},
     numeric, method= bvp[midrich], abserr = 1e-10, output= operator
):

psi:= unapply(eval(x*f(eta)+H(eta), A1), (x,eta)):
plots:-contourplot(psi(x,eta), x= 0..6, eta= 0..6);

interface(rtablesize=16);

I took VV's algorithm for the triangle, rewrote it in a pure matrix form, and made it more efficient. This version can generate 6000 points in 16 milliseconds. Can Mathematica beat that?

restart:
     
TypeTools:-AddType(Point, [realcons$2]);
TypeTools:-AddType(Triangle, [Point$3]);

SampleTriangle:= proc(T::Triangle, n::posint)
local A:= LinearAlgebra:-RandomMatrix((n,2), generator= 0..1., datatype= float[8]);
     evalhf(
          proc(A, n)
          local k;
               for k to n do
                    if A[k,1]+A[k,2] > 1 then
                         A[k,1]:= 1-A[k,1];
                         A[k,2]:= 1-A[k,2]
                    end if
               end do
          end proc
          (A,n)
     );
     <<[1$n]> | A>.< <T[1]> | <T[2]-T[1]> | <T[3]-T[1]>>^+
end proc:

I generate the same triangle as VV:

T:= CodeTools:-Usage(SampleTriangle([[-1,0],[0,4],[2,1]], 6000)):
memory used=0.85MiB, alloc change=0 bytes, cpu time=16.00ms, real time=26.00ms
plot(T, symbol= point, style= point, scaling= constrained);

The plot is the same as VV's except that it has 6000 points instead of 4000.

Here's an example of how to continue a main process after only one of its threads has finished a computation. This uses the basic Threads model rather than the Threads:-Task model.

I'll assume that each of your procedures g, h, i has some sort of main loop so that it can periodically check a boolean variable as a while conditon. The boolean variable is global to the procedures and is called continue, and it's initially set to true. The first thread that finishes sets it to false and returns the correctly computed value. This signals the other threads to return NULL.

In my example, the procedures g, h, i are actually the same procedure. But there would be no essential difference if they were different.

M:= module()
local
     continue:= true,
     ghi:= proc(n::posint, id::integer)
     local k, s:= 0;
          for k to n while continue do
               s:= s+k
          end do;
          if not continue then
               userinfo(1, M, "Thread", id, "aborting; k =", k, ".");
               return
          end if;
          continue:= false;
          userinfo(1, M, "Thread", id, "finished first.");
          s
     end proc,
     ModuleApply:= proc()
     local k,f;
          continue:= true;
          Threads:-Wait(seq(Threads:-Create(ghi(10^5, k), f[k]), k= 1..3));
          {entries(f, 'nolist')}[]    #returns all (very likely just one) non-NULL return values
     end proc
;
end module:          

infolevel[M]:= 1:
M();

ghi: Thread 2 finished first.
ghi: Thread 3 aborting; k = 30316 .
ghi: Thread 1 aborting; k = 47450 .

5000050000

If you need me to show you explicitly how to do this with three different procedures, let me know.

Note for the Threads experts: There is a small chance that more than one procedure will finish at the same time. Nothing bad happens in this case. The procedures will each claim to have finished first, and all return values will be returned.

Here's a naive implementation of a sampler that generates a uniform distribution.

Download SampleTriangle.mw

The specification of options for method= lsode follows an arcane protocol different than that of dsolve's other methods: Some of the options, including mxstep, must be passed in an array (not an Array). The size of the array must be (at least) 21 plus the number of ODEs after the system is converted to a first-order system. The position in the array of the options is documented at ?dsolve,numeric,lsode,advanced. Here's an example:

ode:= diff(y(x),x$4)=y(x):
ic:= y(0)=1, D(y)(0)=2, (D@@2)(y)(0)=3, (D@@3)(y)(0)=4:

#The size of the array must be 21 + #ODEs after conversion to a 1st-order system.
C:= convert(Array(1..25), array):    #The conversion to array is necessary.
C[2]:= 1e-5:   #Absolute error
C[10]:= 10000:  #mxstep: max # of steps allowed

Sol:= dsolve({ode,ic}, numeric, method= lsode, ctrl= C):

Your third question was

when I use lsode which way of numerical solution is applied (Euler,midpoint, rk3, rk4, rkf, heun, ... )?

Hmm, if you don't already know the answer to that, then why did you choose to use lsode---the most complex and difficult-to-use of dsolve's numerous methods? The short answer can be found at ?dsolve,numeric,lsode where it says "an implicit Adams method of evaluation [using] functional iteration (no Jacobian matrix is involved)."

The following procedure generalizes the Answers by Kitonum and One Man so that they'll work for any expression:

Display:= proc(e::algebraic)
     local r:= eval(e);
     print(nprintf("%a = %Q", e, subsindets(e, name, n-> nprintf("%Q", eval(n)))) = r);
     r;
end proc:

To use it, the desired expression must be entered with single quotes around every variable:

(a,b):= (3,4):
c:= Display('a'+'b'):

Now, the value of c has been assigned 7.

Of course, if you use this a lot, you may want to shorten the name Display to a single character.

The code is simple enough to analyze. The relevant code can be seen by

kernelopts(opaquemodules= false):
interface(verboseproc= 2):
eval(combinat:-permute1);

Note that I don't display the actual code here out of respect for copyright issues.

The code is similar to your procedure and only slightly more complicated. I believe that the significant difference is that the Maple code uses a remember table. That means it requires memory for a list of all the permutations, plus lists of all the permutations of all "tails" of the original input list.

Instead of break, put a readstat statement in your code, like this:

for k to 10 do
     if k=5 then readstat("Press ENTER or RETURN to continue:") end if
end do;

For more flexibility, you can start the debugger, like this:

Stop:= ()-> [][]:   #dummy procedure for stopat
stopat(Stop):   #Starts debugger when Stop is executed.

for k to 10 do
     if k=5 then Stop() end if
end do;

If you use sum, then L(i) is evaluated before i is given numeric values. This is called premature evaluation. There are many solutions. The simplest and most efficient is to use add instead of sum. With the special sequencing operators add, mul, and seq, the first argument is not evaluated until after the index variable is replaced by the appropriate values, and it is re-evaluated for each such value.

Premature evaluation is not a bug of sum. For the vast majority of Maple commands, including sum, arguments are fully evaluated before being passed. The evaluation is only condidered premature if it is done before the user wants it to be.

Whatever result you were trying to copy into your post didn't copy. But assuming that it was of the form Sol:= RealRange(a,b), then the b can be extracted as op(2,Sol). If b is of the form Open(c), then c can be extracted as op(op(2,Sol)) or op([2,1], Sol).

ListTools:-PartialSums(A);

It's not necessary that the message be less than N. If that were true, it would severly limit the size of messages which would be feasible to encrypt. Rather, you divide the message into chunks such that each chunk is less than N, and then you encrypt the chunks. It's convenient to use a fixed chunk size of an integer number of bytes. Using your generation technique, I will compute an appropriate chunk size:

(u,v):= 'nextprime(RandomTools:-Generate(integer(range= 10^100..10^(110-1))))' $ 2:
N:= u*v:
C:= iquo(ilog2(N), 8);

It follows that any chunk of 89 bytes or less can be represented as an integer less than N.

I can't test it without the code for RECTpred, RECTpredOR, RECTpredOL. But why aren't you using Threads:-Seq? That's much easier to use and manage than Threads:-Create.

The library commands abs and ceil are highly symbolic, and they're not thread safe (see ?index,threadsafe). They can't be used in code that uses Threads. If their arguments are strictly numeric, then you don't need the symbolic aspects, and you can write trivial thread-safe versions of them:

Abs:= x-> `if`(x < 0, -x, x);

Ceil:= proc(x)
local t:= trunc(x);
     `if`(x=t, x, `if`(x<0, t, t+1))
end proc: