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.

The following worksheet solves your original problem using Statistics:-NonlinearFit. Then it solves the same problem using NonlinearFit with procedural input. You can rewrite the procedure for your more-complicated model.

Download NonlinearFit_proc_input.mw

The default symbol for the imaginary unit in Maple is uppercase I, not the lowercase i that is more usual in mathematics and which you are using above. That default can be changed to any other symbol:

interface(imaginaryunit= i);

But beware that whatever symbol you use, you won't be able to use it as a variable. The use of i as a variable is so common in most programming languages that Maple decided to use the less commomly used symbol I for the imaginary unit.

Form-conversion commands often need to be composed. The evalc command will take most expressions and express them as the sum of their real and imaginary parts. That is, the expression is put into cartesian or rectangular form. When polar is passed a single expression in this form, it converts it to polar form. Then evalf converts to decimal form.

polar(evalc(z));

     polar(sqrt(2), Pi/4)

evalf[4](%);

     polar(1.414, 0.7854)

The command randpoly might be a place to start:

randpoly(x);
randpoly(x, degree= 2);

etc.

Use literally the word undefined. Maple knows what this means.

A simpler way to solve the problem is to pass solve an inequality that forces the solution to be positive. This avoids the need for select, is, or another command line.

f:= gamma*sqrt(4) * sqrt(omega^2*C2^2*R4^2 /
     (C2^4*R4^4*beta^2*gamma^2*omega^4+C2^2*(1+gamma^2*(beta+1)^2-2*gamma) *
      R4^2*omega^2+1)):

solve({diff(f,omega), omega > 0}, omega) assuming positive;

     {omega = 1/(C2*R4*sqrt(beta*gamma))}

I have a hunch that the above is more robust than the flaky and ancient is because it uses the more modern and more sophisticated package SolveTools:-Inequality. However, I also know that it's still easy to confuse solve with inequalities, especially when square roots are involved.

An equivalent command is

solve(diff(f,omega), omega, useassumptions) assuming positive;

It's important for you to understand that the roots that you originally computed with RootFinding:-Isolate (or, equivalently, fsolve) are accurate to 10 digits. The large errors are introduced by your attempt to compute the residuals. Kitonum's Answer may give you the false impression that you need to increase Digits to compute the roots, although that wasn't his intention. No, you only need to increase Digits for the computation of the roots if you also want to compute the residuals, which isn't a necessary part of the computation of the roots.

To say that a matrix is "augmented", in this context, means that the right-side vector b has been appended as the rightmost column of the matrix. To achieve this, you simply need to make the command

IterativeApproximate(<A|b>, ...);

So, the b is incorporated into the first argument. In this case, don't use b as the second argument.

I get this solution:

     Vector[column](4, [3.98912386882810, 2.99292198093478, 2.99317583909596, 2.99337057407024])

which has residual

LinearAlgebra:-Norm(A.% - b);

     0.103910640752635786e-2

which is very close to your specified tolerance of 1e-3.

The above is a response to your originally posted Question, using b = <1,1,1,1>.

Responding to your reformulated Question, using b = <1,0,0,0>: The solution that you got seems correct to me. The residual is 0.235228835640644007e-3, which is well less than your specified tolerance.

It's important for a beginner to note that print doesn't "give" anything that can be stored for later use. It's often difficult for a beginner to understand this. The print command only displays something on the screen; it's rarely used in formal Maple programming.

Here's how to create A in one line of code and then create C in a second line of code:

Download Matrices.mw
A:= Matrix(
     [[0$7], [0$6], [a,-d,a,-d,a], [-a,d,d,a], -[d$3], [d,d]],
     scan= band[0,5],
     shape= antisymmetric
);

C:= [seq(seq(eval(A, {a= i, d= j}), i= [-1,1]), j= [-1,1])];

Now C is a list of four matrices. The matrices can be accessed as C[1], C[2], etc. Note that for loops can often be replaced by seq expressions, and it's usually beneficial to do so.

Example:

A:= Array([4,3,7]);