For more reliabled performance, it is best to not create formulas first and use := to assign values to the variables.  Either assign the variables first and then combine them in formulas, or use eval() to specialize a symbolic formula.

Either

z := 10*Unit('km');
s := 12*Unit('cm');
lambda := 600*10^(-9)*Unit('m');
theta := s/z;
Diam := 1.22*lambda/theta;        

or

theta := s/z;
Diam := 1.22*lambda/theta;
eval(Diam, {lambda = 600*10^(-9)*Unit('m'), s = 12*Unit('cm'), z = 10*Unit('km')});
 

Doing it in the other other, like in your worksheet, works enough of the time that you end up thinking it's okay, but things get evaluated differently and it ends up being problematic when working with more complicated elements like units

P.S. as of Maple 2018, you can just use with(Units): and it will automatically load the best units handling environment.

There haven't been any significant changes made to the top level isprime command in a long time and not even small ones from Maple 18 to 2017. After the initial trial divisions in isprime, Maple does indeed call mpz_probab_prime_p for gmp_isprime, it just chooses its own value for the 'reps' argument.

These particular forms of RootOf is used internally by a few polynomial solvers.  They aren't recognized by everything that recognizes RootOfs and so they aren't emitted in output by any documented commands outside of a few things in RegularChains.

In general, as a user-level command, RootOf is intended for constructing algebraic numbers for use in doing computation in algebraic number fields (see: https://www.maplesoft.com/support/help/maple/view.aspx?path=evala) and not for solving equations.  With few exceptions, it's always better to call the top level solvers (solve, fsolve) for that.  The code in the question is basically a convoluted way of doing:

If you put your system in matrix form, you can call LinearAlgebra:-LinearSolve with the option free=S to have the free variables be called S[1], S[2], etc.  You can use LinearAlgebra:-GenerateMatrix to put your equations into matrix form if you don't have them there already.  Of course, this depends heavily on what your system of equations looks like.  If it's very sparse or has a lot of parameters, you're probably better off trying to post process the output of solve where the free variables will all show up as xn = xn somewhere in the solutions.

Because of all the parameters, just the Row Eschelon Form of this matrix takes 2 gigabytes to store.  That would put the solution at maybe 30 gigabytes.  Due to the large number of parameters, this problem asked for a very different approach from just trying to write down the closed form of the solutions.

You want to format a string, and then use that string to fill the TextArea

y3 := sprintf("%-2.5s:        %2.5s:     %2.5s\n  %-2.5g:   %-2.5g       %-3.2g\n", z, y, x, y1, x1, z1);
DocumentTools:-SetProperty(TextArea0, value, y3);

This is arbitrary and documented in the sqrt help page:

Try right clicking on the equation choosing "Copy As" and selecting "MathML".  Newer versions of Word knows how to eat MathML.

c := ColorTools:-ColorsFromImage("apples.jpg", number = 15)

produces a managable number of colors -- you can easily filter out the whites and greys by hand.
Or automatically:

ColorTools:-Swatches(select((s) -> :-Luma(s) < .6 , c))

You get most of the way there by applying loglogplot then dualaxisplot.

p := plots:-loglogplot(x^3, x = 0 .. 2);
plots:-dualaxisplot(p, p);

When your system includes inequalities, solve implicitly assumes the variables in those inequalities are real.  If you want to be as careful as possible solving polynomial systems of inequalities, you should try the command: SolveTools:-SemiAlgebraic directly, though solve may be calling it in your case (I cannot tell for sure because you have posted a picture of your problem rather than commands that can be cut&pasted).

In general, if f, g, and h are polynomials, you can use the command SolveTools:-SemiAlgebraic to get a mostly out of the box solution.

For this problem, you might get better results by forming the descriminant for the quadratic polynomial in x and solve for the various cases.

fsolve only works on "square" problems where the number of equations matches the number of variables.

If you want to solve a problem with a symbolic parameter (V), you will need to use solve.  Solve works fine on systems with floating point numbers in them, but you will be able to control your results a little better, if you convert them to exact rational numbers first:

eqs := convert([eq1, eq2, eq3, eq4, eq5, eq6], rational);
vars := {c1, c2, c3, phi1, phi2, phi3};

sols := [ solve(eqs, vars) ];

I get 8 solutions.

Try removing some uneval ticks ', there are different evaluation rules inside procedures.

Unfortunately for your application, the numpoints option of implicitplot doesn't actually generate 1000 points on the curve.  The algorithms prune a lot and only generate maybe 100 points for display.  Even if they did, you wouldn't get the same x values for each y so they wouldn't line up in a neat a by x grid anyway.

Here is a way to get all of the x-y points in a Matrix for each value of a:

R0:=gamma+2*ln(2)+Re(Psi(1/2+(1+I)*tanh((1+I)*x)/((tau+0.25e-1*a)*y)))+ln(y);
tau := 9.975;

for a from 1 to 50 do
    P:= plots:-implicitplot(R0, x = 0 .. 5, y = 0 .. 1,numpoints=1000);
    datatmp := [plottools:-getdata(P)];
    data[a] := < ( map2(op,3,datatmp))[] >;
end do;