@Alejandro Jakubi 

In this case I'm inclined to say that the documentation is wrong. But perhaps the behavior of op is wrong. I definitely see nothing wrong with the bahavior of subs.

@Rouben Rostamian  

Preben's explanation in terms of op is not quite correct, but I didn't say anything before because I didn't think an example like yours would come up. The command op does often probe two levels deep into the expression tree, thereby skipping over accessible subunits. To see the complete tree, use ToInert:

P:= 5*a*b:  ToInert(P);

     _Inert_PROD(_Inert_PROD(_Inert_NAME("a"), _Inert_NAME("b")), _Inert_INTPOS(5))

From this, we see that a*b is a distinct subunit of 5*a*b, hence it is accessible by subs even if it's not accessible by op.

I am going to disagree with Alejandro on this point: I think that the "grey box method" ToInert tells you everything that you need to know to understand this situation, although that stuff about assumptions is quite interesting.

@Markiyan Hirnyk 

I prefer

plots:-implicitplot(
     r^2 = 5 - 4*cos(theta), r= -4..4, theta= 0..2*Pi,
     coords= polar, axiscoordinates= polar, gridrefine= 2,
     thickness= 3
);

@baustamm1 

The whole expression is an operand of itself, as can be seen by using unevaluation quotes:

 op('P');

                                      a*b

Your equations are far, far too complicated to solve symbolically. If you supplied numeric values for t and b1, you may be able to solve numerically for gam0, gam1.

By the way, it's a very bad coding practice to use t both as a free variable and a variable of integration.

@jonlg Yes, I realize now that I was overly hasty in asking you for the parameters, because they were in the first attached worksheet. My first reaction was just to look in the code in the main post.

I am working on a solution for you. It's a strange problem that I don't understand, but I think that I have a way to workaround it even if I don't understand it.

@acer 

Your `print/%mixed` is not working for me. Using lprint, I see that %mixed(17/6) is generating the correct function call, %mixed(2, 5/6). But the prettyprint is just giving me NULL. I'm using Maple 16.02, Standard, on Windows 7.1 Home Premium; no initialization file. My Tools->Options->Display settings are Input=Maple Notation, Output=2D Math, Tysetting Level = Maple Standard.

@Kitonum +1 for the excellent 3d cutaway plot.

@dstaple You wrote:

The ideal solution would be equivalent to fsolve(...) or evalf(solve(...)) in the case where the number of equations matched the number of unknowns.

So, are you saying that if the solutions from DirectSearch:-SolveEquations were reformatted so that they were in the same format as fsolve solutions, then that would be an acceptable solution to you? If so, a procedure to do that reformatting would be trivial to write.

@dstaple 

Okay, I understand now that you are looking for a finite set of system-satisfying purely numeric assignments to the variables.

You wrote:

It's possible that this would find a solution in some situations where fsolve would fail; it depends if fsolve already tries this technique as one of its internal algorithms, or if it goes immediately for numerical methods.

As the error message that you got indicates, fsolve will always fail when the number of equations differs from the number of unknowns.

@Markiyan Hirnyk 

Of course, a symbolic solution to a system of nonlinear equations is rarely possible. But I think that you are criticizing my solution because it contains complicated RootOf expressions. A solution expressed in terms of RootOf is still a useful solution. It can be numerically evaluated with evalf, and hence it can be plotted. Here is the most complicated RootOf expression from the original system, the solution for y1 in terms of x2, plotted.

S:= [seq(eliminate({eq||(1..3)}, V), V= combinat:-choose({x1,x2,y1,y2}, 3))]:
Y1:= op([3,1,2,2], S);
plot(Y1, x2= 1..5);

Anyway, in this case you can ignore the RootOf solutions and still have a complete set of solutions with one equation relating every pair of variables.

Please post your complete code. I cannot run your code without numeric values for k, l, m, and r.

@tira 

I corrected all your syntax errors---all the ones mentioned by Preben and several others. Also, I made all the integrals numeric. Then I tested a single point, like this:

evalf(u(1., .4));

(You should always check a single point before plotting a complicated function. If there are any free variables left in the answer, then you know that there are still syntax errors.) The result that I got was

Float(undefined) + 0.*I.

(I was surprised by how quick the answer was returned.) This may indicate that at least one of your integrals does not converge. So I changed the starting point for all the integrals from 0 to 1 to avoid singularities. I still got the above answer. So, I guess that there are other singularities.

I would post the corrected worksheet, but I lost it when I stepped on the power switch of my computer, and I don't have time right now to redo it.

Your delta should be Dirac.

What does f[z1] mean? Is f the same as F?

@tira 

I am not an expert on pdsolve(..., numeric). I do know that it is often wrong. However, whenever I've found it wrong, there have been telltale signs. I do not see those signs in this case. If the situation is that changing a single sign causes you to get the solution that you actually expect (from a done-by-hand symbolic solution), then I strongly suspect that there was a sign error made in that done-by-hand solution.