map(collect, M, [a1, a2]);

To the best of my ability to interpret what you typed, you meant what I have below. But I am still unsure what you meant by ^() (an empty exponent). Note that square brackets cannot be substituted for parentheses in Maple. If I have transcribed your equation correctly, then it can be easily solved with fsolve.

Download fsolve_int.mw

At the end of the computation in Part1, include the commands

try FileTools:-Remove("MyFile.m") catch: end try;
save var1, var2, var3, "MyFile.m";

where var1, var2, var3 are the actual variable names that you want to pass, and there can be any number of them. Make sure that the filename ends with the characters .m. Start Part2 with the command

read "MyFile.m";

Hmm, am I missing something here? Why not just pick values for A and B and plot it?

X:= eval(A*exp(-a*t)+B*exp(-b*t), [a= 1, b= 2, B= 1]):
plot(
     [seq(eval(X, A= k), k= seq(4^k, k= -2..2))], t= -2..2,
     legend= [seq(A= 4^k, k= -2..2)]
);

Addressing your aside only, there is in Maple 18 the InertForm package.

InertForm:-Display(InertForm:-Parse("30 = 6*5"));

The output of this does not display correctly in MaplePrimes. In a Maple 18 worksheet, it looks just like 30 = 6*5. A more cryptic but programatically more accesible wasy to use this is

InertForm:-Display(30 = `%*`(6,5));

Computer algebraically speaking, every infix operator has a prefix form obtained by enclosing it in backquotes. If the prefix form operator is preceded by a percent sign inside the backquotes, that makes it inert.

Okay, here it is coded just with basic for loops:

MyAllNonZero:= proc(M::Matrix)
local m,n,i,j;
     (m,n):= op(1,M); #Get Matrix dimensions.
     for i to m do
          for j to n do
               if M[i,j] = 0 then return false end if
          end do
     end do;
     true
end proc:
         
recu:= proc(C::list(Matrix))
local k;
     for k to nops(C) do
          if MyAllNonZero(C[k]) then return k end if
     end do;
     0
end proc:

It's straightforward list processing. I don't see recursion helping.

recu:= proc(C::list(Matrix))
local k;
     for k to nops(C) do
          if ArrayTools:-AllNonZero(C[k]) then return k end if
     end do;
     0
end proc:

C:= [seq(LinearAlgebra:-RandomMatrix(2,2, density= .5), k= 1..1000)]:
recu(C);

     10

C[1..10];

recu(C[1..9]);

     0

I challenge all other entrants to beat this for time efficiency:

CartProd:= proc(L::list(list))
local S, _i, V:= _i||(1..nops(L));
     [eval(subs(S= seq, foldl(S, [V], (V=~ L)[])))]
end proc:

(m,n,p):= (2,3,3):
map[3](Matrix, m, n, CartProd([[$0..p-1] $ m*n]));

Iterators are memory efficient if the created objects are used and then immediately destroyed. It defeats the purpose if you need to hold all the created objects in memory at once.

Keeping It Simple...

You can use simplify with side relations, which is often more powerful than algsubs.

simplify(A, {Y});

In the future, please copy just your input to your posts, so that a reader can cut-and-paste a single block of code.

Change plots:-pointplot to plots[pointplot]. You must be using a very old version of Maple---from before plots was converted to a module.

It works! Using the simplify with side relations as I described a few hours ago produces the form that you are looking for. To avoid confusion in the future, do not refer to this form as a factorization of the polynomial. I don't know what to properly call it, but it's definitely not a factorization.

Using T2 from your worksheet divide.mw, I did

#For every possible pair of variables, define a new variable as their difference:
Subs:= {seq(`-`(P[])= cat(P[]), P= combinat:-choose({p||(1..4)}, 2))};

simplify(T2, Subs);

#Revert to the original variables:
subs((rhs=lhs)~(Subs), %);

The same process also worked with the final polynomial in your worksheet cubic4.mw, producing a much longer result. Remember to change p||(1..4) to p||(1..6) in the first line (the definition of Subs).

Maple can communicate via sockets. See ?Sockets .

Use $include instead of read. Note that the $ must be in the first column. See ?$include . The problem with read is that the code is read and executed.

You simply need to give a domain restriction to fsolve.

fsolve({F,Fw,Fk,Ft}, {w,k,ki,T} =~ 0..infinity);

Here's a recursive procedure for it, using option remember, though I wouldn't say that it's any more efficient than your loop.

inc1:= ex-> subsindets[flat](ex, indexed, x-> op(0,x)[op(1,x)+1]):
mytest:= proc(C::posint)
local k;
option remember;
     (inc1(thisproc(C-1)) + w[2]*mul(s[k], k= 3..C))*s[1]
end proc:
mytest(2):= s[1]*w[2]+w[1]*s[2]:
mytest(3):= inc1(mytest(2))*s[1]:
mytest(6);