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

If your objective function is a procedure, then the constraints must also be procedures. Your current constraints are expressions. So your command should be like this:

sol:= Optimization:-NLPSolve(
     f, {}, {(x1,x2,x3,x4)-> x1+x2+x3+x4-65}, 0..20, 0..15, 0..20, 0..15, minimize
);

The Units subpackages (Standard in this case) overload the arithmetic operators to behave like this. You can temporarily revert any operator to its default meaning by giving it the :- prefix and using it in prefix form. So, for this case, you could do

eq:= :-`+`(2.*Unit(m), d) = 5.*Unit(m);

     (the equation will display normally: the + will be infix.)

Now you can pass this equation to solve:

solve(eq, d);

e:= w^2/(w^4+2*w^2+1):
expand(numer(e)/w^2)/expand(denom(e)/w^2);

The first expand isn't really needed in this case.

Use 1. * Unit(m). It's not really meaningful to have an exact integer number of something with a unit.

The Hermite form of a matrix of polynomials is a good way to analyze the rank and how the rank depends on the values of the variables.

The following worksheet will not display on MaplePrimes, but please download it.

HermiteForm.mw

Here is the worksheet without output:

Download HermiteForm_no_output.mw