Assuming that you have no interest in the generated file being human readable, then the easiest thing to do is to write and read the matrices in ".m" format. This is Maple's internal format. Essentially the only difference from what you have above is that the format code is %m rather than %a. This format can be used for any Maple data structure, not just matrices.

Download percent-m_format.mw

This is like Kitonum's procedure, but it only displays non-isomorphic graphs, except for n=4. For that n, there is a bug in GraphTheory:-IsIsomorphic. So, for n=4, the procedure displays all 64 "raw" graphs. This procedure takes about 4 minutes for n=6. I wouldn't recommend running it for n=7 or greater. Note that there are over 2 million raw graphs for n=7.

AllGraphs:= proc(n::posint)
uses C= combinat, GT= GraphTheory, LT= ListTools, P= plots;
local  
     Gs:= map2(GT:-Graph, n, [C:-powerset(C:-choose({$1..n},2))[]]),
     NP:= plot(()-> 0, -1..1, color= white, axes= none),
     ng, ns
;    
     if n <> 4 then
          Gs:= map(`?[]`, [LT:-Categorize(GT:-IsIsomorphic, Gs)], [1])
     end if;
     ng:= nops(Gs);
     ns:= ceil(evalf(sqrt(ng)));
     P:-display(Matrix(ns, [GT:-DrawGraph~(Gs, style= circle)[], NP $ ns^2-ng]))     
end proc;

Is your goal to re-import the Matrix into Maple? If so, then use the save command. For example:

A:= LinearAlgebra:-RandomMatrix(31);
save A, "C:/users/owner/desktop/MyMatrix.txt";

To get the Matrix back into Maple, perhaps in another session, do

read "C:/users/owner/desktop/MyMatrix.txt";

This process can be applied to any data structure, not just Matrices.

If your goal is to create a text file that can be read by another program, use the command ExportMatrix. For example,

ExportMatrix("C:/users/owner/desktop/MyMatrix.csv", A, target= csv);

will export your Matrix as comma-separated values (csv). See ?ExportMatrix for the other formats supported.

If you can get the data in a Matrix, it's then trivial to extract the columns from that Matrix. Let's say that the Matrix is named M. If you want the third column of M, you could do

X:= M[.., 3];

See ?printlevel. In particular, setting

printlevel:= 2:

will let you see the results of statements within two nested loops or within a loop and a conditional statement.

Edit: I changed to value from 3 to 2.

The command is evalf[2], for example, evalf[2](2.5999).

You need multiplication operators for xy, xz, and yz. They should be x*y, etc.

It's better to use subscripted variables x[k] rather than concatenated variables xk when there are an indefinite number of them.

fx:= x-> theta*exp(-theta*x);
prod:= product(fx(x[k]), k= 1..n);

In the command I show below, the part Typesetting:-ms(...) is undocumented, and thus it could change in a future Maple release.

plot(
     x^2, x= -2..2,
     title= Typesetting:-ms("My favorite plot", color= green),
     titlefont= [TIMES,ROMAN,24]
);

(FC,Fr):= selectremove(has, F, C):
(FL,Fx1):= selectremove(has, Fr, L):
for v in [C, L, x1] do k||v:= collect(F||v/v, v) end do;

Verify correctness:

expand(kC*C + kL*L + kx1*x1 - F);

Your function f, as presented above, is garbage to Maple because it's missing a multiplication sign between y^3) and (5+. If you correct that, then the command

F:= [op(f)];

will split f into its seven factors, which can be individually accessed as F[1], ..., F[7].

The command

(f1,f2):= selectremove(type, f, `+`);

will split f into two parts, the first being the factored product of all the parts in brackets, and the second being the factored product of the other parts.

Here's how the length of 9 for x+2*y is counted:

Running dismantle(x+2*y) shows:

SUM(5)
   NAME(4): x
   INTPOS(2): 1
   NAME(4): y
   INTPOS(2): 2

(The numbers that appear in parentheses after the keywords are irrelevant to this discussion.) This shows that the expression has four operands: x, 1, y, 2. The 1 is the implied coefficient of x. Each of those operands has length 1, as you can verify. So, the expression x+2*y is stored in memory as

DA:= disassemble(addressof(x+2*y));

DA:= 16, 18446744562576677534, 18446744073709551617, 18446744562576677566, 18446744073709551619

The 16 is a tag that indicates that this is a SUM:

kernelopts(dagtag= 16);

SUM

The four long numbers are the memory addresses of x, 1, y, and 2. This can be verified by

seq(pointto(k), k= DA[2..]);

x, 1, y, 2

These addesses are, of course, session dependent. So, the expression x+2*y is stored as five words (that includes one word for the dag-tag 16). The recursive length rule, which you quoted, says that this is added to the lengths of the operands, which is 4 (each operand has length 1, as noted above). Adding 5 and 4 makes the length 9.

Note that multiplication by a numeric coefficient is not considered as a separate multiplication operation.

If you simply want the length in characters of an expression e, this can be easily obtained as 

length(convert(e, string));

(as already mentioned by Axel).

For a command that uses algebra rather than nitty-gritty internal representations to measure expressions, look at SolveTools:-Complexity. Another option is to convert your expression to a procedure (just wrap it with unapply) and use SoftwareMetrics:-HalsteadMetrics. This'll measure the complexity in several ways that attempt to account for the psychology of comphrehending expressions. Two rather simplistic metrics are the MmaTranslator:-Mma:-LeafCount mentioned by Axel, and a simple operation count done by codegen:-cost. 

To my mind, common subexpressions shouldn't count towards the complexity of an expression. Thus, my favorite way to measure the complexity of an expression is to pass it to codegen:-optimize with the tryhard option and to then measure the length (in characters!) of the resulting procedure, after compression. This can be done in one line with

length(sprintf("%m", codegen:-optimize(unapply(e), tryhard)));

where e is the expression. For trivial expressions, this returns a result that's a bit high because there's the overhead of the procedure header. But, you're not dealing with trivial expressions. You can normalize for this effect by simply subtracting 31 from the result, which is the count when the raw metric is applied to the simplest possible algebraic expression, 0.

All inequality constraints should go in the first set; in your case, that'd be the empty set. All equality constraints should go in the second set. So the NLPSolve command should be

sol := Optimization:-NLPSolve(
     f, {}, {p1, p2, p3}, 0 .. 1, 0 .. 1, 0 .. 1, 0 .. 1, 0 .. 1, 0 .. 1, initialpoint = [.5, .5, .5, .5, .5, .5]
);

which can be shortened to

sol:= Optimization:-NLPSolve(f, {}, {p||(1..3)}, (0..1)$6, initialpoint= [.5$6]);

This generates a new error message:

Error, (in Optimization:-NLPSolve) matrix dimensions don't match: 3 x 6 vs 3 x 2

This seems to be a bug in NLPSolve. Running trace(f) shows that your procedure f was never called, so the dimensions problem is not in your code.

Use command march (short for Maple archive). Use march(create, ...) to create an archive. Use march(add, ...) to add the .m file to the archive. Use march(list, ...) to list the contents of the archive.

You use a global variable t symbolically; it's never given a value. Thus, the value returned from f is an expression depending on t[1], t[2], and t[3]. The value returned by f needs to be numeric if it's to be used by NLPSolve. If you do trace(f) before calling NLPSolve, you can see this.

Don't set warnlevel to 0 before your code is debugged. That's just asking for trouble. If the warning tells you that the variables are local, then declare them local.

You are using some very old coding styles. Perhaps you are emulating a textbook. It's giving you bad coding habits. There's no need for with(linalg). Change array to Array.

Don't rely on the with command to allow the use of so-called "short form" names inside a procedure. Instead, use a uses clause in the procedure header, or just use the "long form" name.