I wrote a procedure to compute an approximation to the limit that you gave in your most-recent Reply:

lambda(x0) = limit(sum(ln(abs(f'(x[i]))), i= 0..n-1)/n, n= infinity),

where x[i] = f(x[i-1]), x[0]= x0.

(I figured this out from this Wikipedia article on Lyapunov exponents.)

For your sample function f:= x-> exp(x^2*(a-x)), the convergence is very slow for most of the x0 that I tried. Does anyone here know any numeric acceleration techniques for a sequence of this type? It's almost a series. Fortunately, everything except the initial differentiation of f can be done under evalhf.

LyapunovExponent:= proc(
     f::procedure,
     x0::complexcons,
     {epsilon::positive:= 1e-7}, #Absolute error tolerance
     {max_iters::posint:= 2^21}   #Infinite loop prevention
)
option `Author: Carl J. Love, 2016-May-10`;
description
"Computes a numeric approximation to "
"lambda(x0) = limit(sum(ln(abs(f'(x[i]))), i= 0..n-1)/n, n= infinity), "
"where x[i] = f(x[i-1]), x[0]= x0. "
"See: https://en.wikipedia.org/wiki/Lyapunov_exponent"
;
local f1:= D(f);
     evalhf(
          proc(f, f1, x0, epsilon, max_iters)
          local
               S:= 0,
               L:= epsilon,
               newL:= 0,
               x:= x0,
               n  #loop index
          ;
               for n to max_iters while abs(newL-L) >= epsilon do
                    S:= S+ln(abs(f1(x)));
                    L:= newL;
                    newL:= S/n;
                    x:= f(x)
               end do;
               if n > max_iters then
                    error
                         "Didn't converge. Abs. err. = %1. "
                         "Increasing epsilon or max_iters might help.",
                         abs(newL-L)
               end if;
               userinfo(1, LyapunovExponent, "Converged; iterations =", n);
               newL
          end proc
          (f, f1, x0, epsilon, max_iters)
     )
end proc:

Here's your function done with a sample a and x0.

infolevel[LyapunovExponent]:= 1:

f:= exp(x^2*(a-x)): x0:= 2:
LyapunovExponent(unapply(eval(f, a= .1), x), x0);
LyapunovExponent: Converged; iterations = 2036767.
                     -0.0552719105995118365

So, would you like some plots using this? There are five dimensions to work with: three real and two imaginary since both a and x0 can be complex. This function exp(x^2*(a-x)) is your baby; I'm at a loss for coming up with good ranges for a and x0.

Use seq to solve your problem:

solve({seq(eqc[j], j= 1..2*M-1, 2)}, {seq(C[j], j= 1..2*M-1, 2)});

By the way, your title is understandable; however, correct English would be How to solve for a variable number of variables? The word number, in this context, is used for quantities that are necessarily nonnegative integers; the word amount is reserved for quantities that aren't necessarily nonnegative integers. The for is required because the phrasal verb solve for has a different meaning than solve. One solves an equation; one solves for a variable in an equation.

Below, I've simplified and modernized your code. I agree with Mac Dude that unassigning delta is useless because it was never assigned a value in the first place.

I thought that it would be useful to you to be able to see the solution set returned by solve, so I added code to print that out. Then maybe you can figure what's happening to your delta.

derivation := proc(A::Array)
local
     i, j, k, m, s, D, delta, n:= op([2,2,2], A), _D:= Matrix((n,n), symbol= D),
     sols:= [solve(
          {seq(
               seq(seq(add(A[k,j,m]*D[k,i]+delta*A[i,k,m]*D[k,j], k= 1..n), m= 1..n), j= 1..n),  i= 1..n)
          },
          indets(_D)
     )]
;
     userinfo(1, ':-sols', NoName, print('sols'= sols));
     seq(eval(_D, s), s= sols)
end proc:

AS1:= Array(((1..2)$3), {(1,1,2)= 1}):
infolevel[sols]:= 1:
derivation(AS1);

     sols = [{D[1, 1] = 0, D[1, 2] = 0, D[2, 1] = D[2, 1], D[2, 2] = D[2, 2]}]
     Matrix(2, 2, [[0, 0], [D[2, 1], D[2, 2]]])

It's different if I change the solved-for variables from indets(_D) to indets(_D) union {delta}. Then there are two solutions and consequently two matrices.

derivation(AS1);

     sols = [{delta = delta, D[1, 1] = 0, D[1, 2] = 0, D[2, 1] = D[2, 1], D[2, 2] = D[2, 2]},
          {delta = -1, D[1, 1] = D[1, 1], D[1, 2] = 0, D[2, 1] = D[2, 1], D[2, 2] = D[2, 2]}]

     Matrix(2, 2, [[0, 0], [D[2, 1], D[2, 2]]]), Matrix(2, 2, [[D[1, 1], 0], [D[2, 1], D[2, 2]]])

If you use add instead of sum, then you can use a step of 2.

X:= m-> (1-cos(m*Pi*x/2/a))/2;
w:= unapply(add(C[i]*X(i), i= 1..M, 2), x);

This assumes that M has a definite numeric value. If you want M to remain indefinite, then let me know. It's fairly easy to re-arrange it for that.

Here's a very large start to the problem. You can also make 2-D plots by changing [x,y,z] in the plot command to [x,y], [x,z], or [y,z].

I don't exactly know what you mean by "code for bifurcation by varying tau vs x(t)." Perhaps you could point me to an article about it.

restart:
ODEs:=
     diff(x(t),t) = y(t)-b*x(t)^3+a*x(t)^2-z(t-tau)+J,
     diff(y(t),t) = c-d*x(t)^2-y(t),
     diff(z(t),t) = r*(s*(x(t)-beta)-z(t))
:
ICs:= x(0)=2.6, y(0)=1.7, z(0)=0.8:
params:= {a= 3, b= 1, c= 1, d= 5, s= 4, beta= 1.6, r= 0.6e-2, J= 3.0}:

#Plotting parameters:
taumax:= 10:
tmax:= 20:
ntau:= 5:

Sol:= tau-> dsolve(eval({ODEs, ICs}, params union {:-tau= tau}), numeric):

plots:-display(
     seq(
          plots:-odeplot(
               Sol(tau), [x,y,z](t), t= 0..tmax,
               numpoints= trunc(50*tmax), color= COLOR(HUE, tau/taumax)
          ),
          tau= 0..taumax, taumax/ntau
     ),
     thickness= 0
);

Excellent Question. Here's a workaround. The explanation for why what you tried didn't work follows.

F_if:= unapply(evalindets(F(x,y), specfunc(anything, piecewise), convert, `if`), (x,y));

     F_if:= (x,y)-> if x < 0 then if y < 0 then 0 else 1 end if else if y < 0 then 1 else 0 end if end if

That's word-for-word what you were expecting!

Explanation: The command convert(..., `if`) is, as far as I can see, undocumented. However, it's natural that you would try it because the documented command convert(..., piecewise) does the inverse conversion. The code for this command is in the procedure `convert/if`, which can be displayed with

showstat(`convert/if`);

(I don't display the actual code here due to copyright issues.) We see that this is a one-line procedure that relies on Maple's "hackware package": the commands pointto, assemble, addressof, and kernelopts(dagtag= ...) (these commands all have help pages). Even without knowing what those do, it's easy to see that it only works on the top level (see footnote *). The design is flawed in that the procedure doesn't apply itself recursively to piecewise substructures. This can be corrected with command evalindets. That'll extract all substructures of a given type and recursively apply a conversion to them all.

*Expert-level observation: We can see from this procedure that the internal structure of an if statement has its arguments in the same order as a piecewise command.

These Mathematica commands for finite-element method with adaptive multigrid (mesh refinement) are light-years ahead of what's possible for the numeric solution of PDEs in Maple. There are no equivalent Maple commands, so no direct translation is possible.

Of course Maple does provide all the low-level programming capability required to program these methods from scratch. Searching the Maple Applications Center for "finite element method PDE" turns up one hit. There may be other third-party resources on the web.

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.