It works for me. You'll need to post an example of it not working for you.

restart;
p:= "acat:= \"a cat\"";
                        p:= "acat:= "a cat""
writeline("junk.mpl", p);
                               15
restart;
p:= readline("junk.mpl");
                        p:= "acat:= "a cat""
parse(p, statement);
                           acat:= "a cat"
                            "a cat"
acat;
                            "a cat"

The two solutions need to be enclosed in square brackets so that they both can be plotted.

c:= x^2+(y-3)^2=25:
cplot:= [solve](c,y):

Notice the square brackets in the above command.

plot(cplot, x= -5..5, scaling= constrained);

How about this? I just made some modifications to your code. For each cylinder I pick two points and the radius at random. The points are the centers of the endcaps. Then I do a tubeplot of the line segment connecting the points.

with(RandomTools):
with(plots):

x:=Generate(list(integer(range=0..20),100)): # random coordinate for x parameter
x1:= Generate(list(integer(range= 0..20), 100)):
y:=Generate(list(integer(range=0..20),100)): # random coordinate for y parameter
y1:= Generate(list(integer(range= 0..20), 100)):
z:=Generate(list(integer(range=0..20),100)): # random coordinate for z paramete
z1:= Generate(list(integer(range= 0..20), 100)):
R:=Generate(list(float(range=1..3),12)):  # random radius parameter

P:=(x,x1,y,y1,z,z1,R)-> [x+(x1-x)*t, y+(y1-y)*t, z+(z1-z)*t, radius= R]: # equation for one cylinder
S:=(x,x1,y,y1,z,z1,R)-> tubeplot(P(x,x1,y,y1,z,z1,R), t= 0..1, tubepoints= 20): # 3D plot in space

display3d({seq}(S(x[k],x1[k],y[k],y1[k],z[k],z1[k],R[k]),k=1..12));

If you replace abs(x) with csgn(x)*x, then f(x) is analytic in a disk in the complex plane about x=1, and the taylor command (in M16) produces the expected result. This is equal to abs(x) for real x.

The fact that Maple (16) does not give a Taylor series expansion for sqrt(abs(x)) about x=1 shows that a bug has been corrected, not that a bug has been introduced.

Your original form [[1,2],[3,4]] is not a Vector, it is a list of lists. To get back to that form

convert(A1, listlist);

Addressing your other question, you simply need to take the transpose of A1 before converting:

convert(A1^+, listlist);

Here's another solution, which uses elementary number theory (in particular the Fundamental Theorem of Arithmetic) to substantially reduce the search space, and thus get the answer in under 0.1 seconds.

restart;
st:= time():

If we multiply each of a,b,c,d by 100, then the sum is 711, and the product is

N:= 711*100^3:

and now the problem is strictly about positive integers. Henceforth, a,b,c,d will only refer to these integerized values.

Each of a,b,c,d must be a divisor of N, so let's look at the prime factorization of N:

ifactor(N);
                           6    2    6     
                      (2)  (3)  (5)  (79)

There is a prime factor that only appears to the first power. That will speed this up a lot. Let's get the divisors of N, which will be our search space.

divs:= numtheory:-divisors(N):

Since the sum is 711, we restrict the search to divisors less than 711.

divs:= select(`\<`, divs, 711)

(Ignore the backslash in the above line. Maple ignores it, and the MaplePrimes editor in my web browser won't display the rest of the line without it.)

Exactly one of a,b,c,d must be a multiple of 79---we'll make it a---so we partition the remaining divisors based on that.

d79, divs:= selectremove(x-> x mod 79 = 0, divs):

Now loop through these sets. After each x is chosen, the remaining solutions must divide N/x, so we reduce the search space accordingly in each inner loop.

proc()
     local Na, divsa, divsb, Nb, a, b, c, d, sol;
     for a in d79 do
          Na:= N/a;
          divsa:= divs intersect numtheory:-divisors(Na);
          for b in divsa do
               Nb:= Na/b;
               divsb:= divsa intersect numtheory:-divisors(Nb);
               for c in divsb do
                   if a+b+c+Nb/c = 711 then sol[a,b,c]:= {a,b,c,Nb/c} fi
               od
          od
     od;
     convert(sol, set)
end proc
();
                          { {120, 125, 150, 316} }

Finally, divide each of those by 100 to get the desired answers. I put each solution in a set so that permutations would be condensed. If there was a solution with repeated values, it would appear as a set with fewer than four members, and we'd have to do a little work to figure out which values were repeated. But, there is no such solution.

time()-st;

                                   0.031

In a separate run, I counted the number of executions of the innermost loop. It was 8520. I also tried the above without separating the multiple-of-79 divisors as special. Then I got 58,938 iterations in 0.109 seconds, so maybe it wasn't worth it to separate them.

A note on coding style: I use an anonymous procedure above, i.e. a procdeure with no name. I invoke the procedure immediately by placing () after its definition.

I changed your example procedure to one where the order of the arguments matters so that the example is more generally applicable.

myproc:= (a,b,c)-> a/b-c:

The `~` command below applies the procedure in "batch" form.

(myproc@op) ~ ([[d,e,f], [g,h,i]]);

Because of operator precedences, all the parentheses in the above command are necessary, unfortunately. Also, the arguments must be in lists [] rather than sets {} because the order matters. I think this command is fairly close to your desired syntax.

The above command is equivalent to the more traditional form

map(myproc@op, [[d,e,f], [g,h,i]]);

The op is used to convert lists, such as [d,e,f], to sequences, such as d, e, f; since myproc expects a sequence of arguments.
                      

Example 1: Factor the first 2⁹ odd 2⁷-bit integers. First, the "old" way:

restart;
N:= 2^(2^7-1)+1:
B:= [seq](N+k, k= 0..2^9-1, 2):
CodeTools:-Usage(map(ifactor, B)):

memory used=1.12GiB, alloc change=212.68MiB, cpu time=18.76s, real time=21.10s

...and now using package Threads:

restart;
N:= 2^(2^7-1)+1:
B:= [seq](N+k, k= 0..2^9-1, 2):
CodeTools:-Usage(Threads:-Map(ifactor, B)):

memory used=1.23GiB, alloc change=0.53GiB, cpu time=33.50s, real time=8.97s

Example 2: Generate 2¹⁴ random 2⁷-bit primes. The old way:

restart;
N:= 2^(2^7-1):
Gen:= RandomTools:-MersenneTwister:-NewGenerator(range= N..2*N):
seeds:= [seq](Gen(), k= 1..2^14):
Primes:= CodeTools:-Usage(map(nextprime, seeds)):

memory used=2.76GiB, alloc change=327.83MiB, cpu time=27.50s, real time=27.50s

Primes[-1];  #Inspect last prime in list...

            197727613027066792031910829881298775459

...and using Threads:

restart;
N:= 2^(2^7-1):
Gen:= RandomTools:-MersenneTwister:-NewGenerator(range= N..2*N):
seeds:= [seq](Gen(), k= 1..2^14):
Primes:= CodeTools:-Usage(Threads:-Map(nextprime, seeds)):

memory used=3.74GiB, alloc change=1.56GiB, cpu time=77.12s, real time=11.00s

Primes[-1];  #...just to verify that the two methods are doing the same computation.

            197727613027066792031910829881298775459

Records are mutable, as the following shows:

restart;
A:= Record(a=1):
B:= Record(a=1):
T[A]:= 1: T[B]:= 2:
eval(T);

   table([B = 2, A = 1])

Compare that with what happens with an immutable structure, in this case an unevaluated function call (named `record`, just for fun):

restart;
A:= record(a=1):
B:= record(a=1):
T[A]:= 1: T[B]:= 2:
eval(T);

   table([record(a = 1) = 2])

But that doesn't mean that they absolutely can't be used with remember tables. If the records being compared are address identical, not just copies which "just happen to be" equal, then the remember table will recognize them. If just a few cases slip through unrecognized because they just happen to be equal, it won't ruin your recursion. On the other hand, if your algorithm relies on making copies of records (not just copying the address, or pointer), the remember table won't work.

Cache tables, remember tables, and just-plain tables all function like the tables in the examples above. The only difference that cache tables have is that the amount of memory that they use can be easily controlled.

restart;
`convert/range`:= proc(ineqs::{set,list})
# Convert a pair of simple inequalities to form x= a..b.
   local ab:= indets(ineqs, numeric);
   indets(ineqs, name)[]= min(ab)..max(ab)
end proc:
n:= 1..3:
s1:= (x-1)*(y-x):
s2:= (7-y)*(1-x):
s3:= (x-y)*(y-7):
feas:= [solve](({s||n} >=~ 0) union {-2<=x,x<=3, 0<=y,y<=11}, {x,y});

      [{y = x, 0 <= x, x <= 1}, {y = 7, 1 <= x, x <= 3},
        {x = 1, 1 <= y, y <= 7}]

z:= add(sqrt(s||i), i= n):
m:= 1..nops(feas):
zr||m:= seq(eval(z, indets(feas[i], `=`)), i= m):
seq(maximize(zr||i, convert(indets(feas[i], `<=`), range)), i= m);

                      (1/2)   (1/2)   (1/2)
                     7     , 8     , 9   

So the answer is sqrt(9), or 3.

In your example, you don't need to create the Array until you know how many entries it will need. Presumably blocks is a list with the same number of entries as you want in the Array. That number can be accessed as nops(blocks). So right after the line where you create blocks, you can put the line

c:= Array(1..nops(blocks));

You can eliminate the for loop by including an initializer in the above statement, and putting it after the definitions for e and m:

c:= Array(1..nops(blocks), i-> fme(convert(blocks[i], bytes), e, m));

and not printing c until after the procedure returns.

Here's my solution:

> restart;
    invtau:= proc(n::posint)
   # Inverse tau: Returns the smallest positive integer r such
   # tau(r) = n, where tau(r) is the number of positive integer
   # divisors of r (see ?numtheory,tau).
       local pair, r:= 1, r_prime:= 1, exponent;
       for pair in ListTools:-Reverse(sort(ifactors(n)[2])) do
           exponent:= pair[1] - 1;
           to pair[2] do
               r_prime:= nextprime(r_prime);
               r:= r*r_prime^exponent
           od
       od;
       r
  end proc:

> st:= time():
> [seq]([n,invtau(n)], n= 1..100);

[[1, 1], [2, 2], [3, 4], [4, 6], [5, 16], [6, 12], [7, 64], [8, 30], [9, 36], [10, 48], [11, 1024], [12, 60], [13, 4096], [14, 192], [15, 144], [16, 210], [17, 65536], [18, 180], [19, 262144], [20, 240], [21, 576], [22, 3072], [23, 4194304], [24, 420], [25, 1296], [26, 12288], [27, 900], [28, 960], [29, 268435456], [30, 720], [31, 1073741824], [32, 2310], [33, 9216], [34, 196608], [35, 5184], [36, 1260], [37, 68719476736], [38, 786432], [39, 36864], [40, 1680], [41, 1099511627776], [42, 2880], [43, 4398046511104], [44, 15360], [45, 3600], [46, 12582912], [47, 70368744177664], [48, 4620], [49, 46656], [50, 6480], [51, 589824], [52, 61440],
[53, 4503599627370496], [54, 6300], [55, 82944], [56, 6720], [57, 2359296], [58, 805306368], [59, 288230376151711744], [60, 5040], [61, 1152921504606846976], [62, 3221225472], [63, 14400], [64, 30030], [65, 331776], [66, 46080], [67, 73786976294838206464], [68, 983040], [69, 37748736], [70, 25920], [71, 1180591620717411303424], [72, 13860], [73, 4722366482869645213696], [74, 206158430208], [75, 32400], [76, 3932160], [77, 746496], [78, 184320], [79, 302231454903657293676544], [80, 18480], [81, 44100], [82, 3298534883328], [83, 4835703278458516698824704], [84, 20160], [85, 5308416], [86, 13194139533312], [87, 2415919104], [88, 107520], [89, 309485009821345068724781056], [90, 25200], [91, 2985984], [92, 62914560], [93, 9663676416], [94, 211106232532992], [95, 21233664], [96, 60060], [97, 79228162514264337593543950336], [98, 233280], [99, 230400], [100, 45360]]

> time()-st;
                             0.031

A:= x^n+y^3.5:
                     
indets(A, anything^float);
                             / 3.5\
                            { y    }
                             \     /
indets(A, anything^symbol);
                                n
                             { x  }
                             
indets(A, `^`);
                             n   3.5
                          { x , y    }
                                   

There is a way to solve this without resorting to commands outside of dsolve. The help page ?dsolve,numeric_bvp,advanced gives advice on how to apply dsolve's options to resolve some error messages given by `dsolve/numeric/bvp`. For messages related to Newton iteration convergence, it suggests the 'continuation' option, which is defined on help page ?dsolve,numeric,bvp (with two examples of its use given on the former help page).

To use the method, one introduces a parameter λ into the problem such that λ can vary continuously from 0 to 1, with λ=0 yielding a relatively easier BVP and λ=1 giving the desired BVP. The parameter can appear anywhere in the problem: in the ODEs, or the boundary conditions, or both. There is some subtlety required in deciding where and how to place λ. I tried this two ways with the problem at hand, and they both worked.

My first idea was to have λ=0 eliminate the convolution of f(x) and g(x), so I applied it to the g(x) term in the first ODE and to the term containing f(x) in the second ODE.

restart;
ode1:= ( (D@@3)(f) + 3*f*(D@@2)(f) - 2*D(f)^2 ) (x) + lambda*g(x) = 0;
ode2:= (D@@2)(g)(x) + lambda*3*10*(f*D(g))(x) = 0;
bcs1:= D(f)(0) = 0, f(0) = 0, D(f)(6) = 0;
bcs2:= g(0) = 1, g(6) = 0;
sys:= {bcs1, bcs2, ode1, ode2}:
dsn:= dsolve(sys, numeric, continuation= lambda);
plots:-odeplot(dsn, [x, g(x)], 0 .. 6, color= black);

I recommend plotting the interval 5.7 .. 6 separately because there is an interesting fluctuation there which is not visible at the g-axis scale of the overall plot.

My second idea was to use λ=0 to "level" the boundary conditions so that everything starts and ends at 0 on the vertical axis. In this case, that simply means changing g(0)=1 to g(0)=λ.

restart;
ode1:= ( (D@@3)(f) + 3*f*(D@@2)(f) - 2*D(f)^2 + g ) (x) = 0;
ode2:= ( (D@@2)(g) + 3*10*f*D(g) ) (x) = 0;
bcs1:= D(f)(0) = 0, f(0) = 0, D(f)(6) = 0;
bcs2:= g(0) = lambda, g(6) = 0;
sys:= {bcs1, bcs2, ode1, ode2}:
dsn:= dsolve(sys, numeric, continuation= lambda);
plots:-odeplot(dsn, [x, g(x)], 0 .. 6, color = black);

Same plot, of course.

You are mixing new and old matrix commands. The new commands all begin with capital letters and tend to be spelled in full rather than abbreviated. Use Matrix and Vector instead of array. Use Adjoint instead of adjoint. Use Determinant instead of Det. See ?Matrix ?Vector ?Determinant etc.