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;

It is often difficult for Maple to do much with symbolic exponents.  If you assign a value to w, you may get some reasonable solutions in terms of the other parameters, but for most values you will probably still get RootOf expressions since the roots of higher degree polynomials aren't generally expressible in other ways.

You want sscanf not sprintf to go from strings:

n := sscanf("123", "%d")[1];

Check out the documentation for dsolve/numeric it should get you started.

It is probably not worth it.  It is very likely to be faster to use expressions directly -- converting to and from functions is just going to add overhead.

There are a couple things that could be slowing this down, but the main one is likely that the Sigma symbol parses into sum() which is much slower than the more straight forward add().  Try defining J this way:

J := unapply(8*Pi^(3/2)*r*R*(add((2*r*R)^(2*i)*pochhammer((1/2)*n, i)
 *pochhammer((1/2)*n+1/2, i)*(add((-1)^j*cos(phi)^(2*j)
 *(add((2*r*cos(phi))^(2.*l)*pochhammer(n+2*i, 2*l)
 *hypergeom([2*j+2*l+1, .5], [2*j+2*l+1.5], -1)*(.5
 *Beta(l+.5, n+2*i+l-.5)-sin(arctan(-Z/sqrt(R^2+r^2)))^(2*l+1)
 *hypergeom([-n-2*i-l+1.5, l+.5], [l+1.5], sin(arctan(-Z/sqrt(R^2+r^2)))^2)/(2*l+1))
/(factorial(2*l)*pochhammer(2*j+2*l+1, .5)*(R^2+r^2)^(n+2*i+l-.5)),
 l = 0 .. 100))/(factorial(i-j)*factorial(j)),
 j = 0 .. i))/factorial(i), i = 0 .. 100)),[n,phi]);

Your worksheet still takes about three minutes to execute, but I killed the other version after 10 minutes.

A couple other answers have the obvious solution, which invovles first creating a new list with the zeros removed.  But the most readble way to do that (imo) has not been mentioned:

remove(`=`, L, 0);

If for some reason you absolutely need to avoid making a new list, you can use ListTools:-FindMinimalElements with a custom comparator that treats 0's as greater than all other elements.

L := [0, 5, 0, 10, 3, 0];
ListTools:-FindMinimalElement(L, (x, y) -> `if`(x = 0, false, `if`(y = 0, true, evalb(x < y))) );

It is probably a pretty rare situation where this is necessary.  Just creating a new list with 0's removed then using min is almost certianly more efficient most of the time.

This is a slightly more generalizable version of Robert's coeffs answer.

vars := {a,b};
eqs := {coeffs((rhs-lhs)(f), indets(f) minus vars)};
solve(eqs, vars);

Assuming your basis polynomials are in P and assuming they are not already a groebner basis, this is what I would do.  Notice that if NormalForm does not return r=0, then it is not possible to write f in terms of P with out a remainder.

> with(Groebner);
> P := [x*y, x^2+y^2]; v := tdeg(x, y);
> G, C := Basis(P, v, output = extended); # C tells you how to write G in terms of P
> psubs := {seq(g || i = add(C[i, j]*p || j, j = 1 .. nops(P)), i = 1 .. nops(G))}:
> f := (x+y)^4; r:=NormalForm(f, G, v, 'Q'); Q; # Q is assigned by side-effect
                             0
                  [   2              2   2   ]
                  [4 x  + 5 x y + 4 y , x , y]
> gb := r + add(Q[i]*g || i, i = 1 .. nops(Q)); # writes f in terms of G
            /   2              2\       2          
            \4 x  + 5 x y + 4 y / g1 + x  g2 + y g3
> pb := collect(eval(gb, psubs), [seq(p || i, i = 1 .. nops(P))]); # and then in terms of P
            /   2              2\      / 2    2\   
            \4 x  + 4 x y + 4 y / p1 + \x  + y / p2