To retain the symbolic constants, use solve instead of fsolve.

There is a way to create true symbolic constants, but you can't get fsolve to express solutions in terms of them; fsolve is strictly numeric.

Two ways. In the first, change g to

g:= proc(t)
local j;
     if not t::numeric then return 'procname'(t) end if;
    ... the rest of g ...

And change the plot command to

plot([g(t)[1], g(t)[2], t= 0..1]);

The second way is to simplify everything:

f0:= t-> (t, 3.9*t*(1-t)):
IFS:= (i,x,y)-> piecewise(i=0, [y,x], i=1, [x,y+1], i=2, [x,y]+~1, i=3, [2-y,1-x])[]/2:
g:= t-> IFS(trunc(4*t), f0(frac(4*t))):
plot([g(t)[1], g(t)[2], t= 0..1]);

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
);