Why have you put m3 in square brackets? That's the problem.

Also with*matrix and with*linearalgebra don't mean anything. Get rid of those.

Try changing Do(HD= HuffmanCoding(HT)) to
Do(HD= [HuffmanCoding(HT)]).

The reason is that I doubt that Do can handle the assignment of an expression sequence with more than one term, which is what the output of HuffmanCoding will be.

Here's an example:

Download DataFrame_string_filter.mw

However, it's a disappointment that the commands select, remove, and selectremove haven't been overloaded to work with DataSeries and DataFrames.

Here's an easy example that you should be able to execute on your own computer.

I'm nearly certain at this point, so I'll make this an Answer:

I think that the last statement of OrderB is of the form print(x). You need to change that to either return(x), return x, or simply x. It doesn't matter which one. A procedure should never use print to return its final value. The purpose of print is to display supplementary information on the screen.

The differential equations and the initial conditions need to be in a single set. So, you need to

dsolve(IVP union ICs, numeric);

Delta:= D[2,1,1] + D[2,1] + D[2,3,3] + D[1,1] + D[3,1,1]:

To use it:

Delta(f)(x,y,z);

Here's a workaround.

Download cosh_integral.mw

Obviously the result is true for negative k also.

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.