@Carl Love 

Mathematically, the flux is a surface integral; it reduces to a double integral. If the integral exists (i.e. the integrand is an integrable function in any sense - Riemann/improper/Lebesgue etc) then there is no math problem. 

In this case it exists. Note that many improper integrals become proper after a change of variables.

The divergence theorem was not used here; and it cannot be used. Actually Div(vv) is not integrable in any ball centered at <0,0,0>.

@acer 

The  "By hand" method is the easiest for a beginner: simply copy the output of Flux(...,inert) and paste it in a new execution group. Then swap  phi = 0 .. Pi    and  theta = 0 .. 2*Pi.
(with no escaped variables problem).

@EugeneKalentev 

That is exactly what evalf does if applied to an inert integral (see ?int/details): the integral is computed numerically even if a closed form exists.

@acer 

It is easy to isolate the integral which Maple cannot compute. It is

f:=sin(phi)^2*sin(r*cos(phi)):
J:=Int(f, phi = 0 .. Pi);

It does not have an elementary primitive (Risch' algorithm is implemented for such functions AFAIK - i.e. not containing algebraic parts).
Your change of variable ==>
IntegrationTools:-Change(J, cos(phi) = t, t);

This an integral of an odd function on a symmetric interval and Maple knows that it is 0.
The simpler change of variable also works (same reason):

IntegrationTools:-Change(J, phi = Pi/2 - t, t);

Notice that cos(phi)=t also works e.g. for phi = 0..Pi/2 (i.e. flux through a hemisphere)

Int(f, phi = 0 .. Pi/2):
IntegrationTools:-Change(%, cos(phi) = t, t);

value(%);

But not for phi = 0..Pi/3 !!

However changing the order of integration works in these cases too, simply because the integral J vanishes (being multiplied by a zero factor).
Of course the method will also fail if theta runs in other intervals e.g. 0..Pi/6.

Best regards,
V.

@ecterrab 

Thank you again.

But I want to mention that the LambertW solution was produced even if the library is not updated (yet).

@ecterrab 

Thank you for the competent analysis.

Could you please also answer to the following questions.

1. Why is dsolve giving in this case different answers? For example, after loading your worksheet, dsolve gave the LambertW solution, not the RootOf one.
2. Is the library contained in Physics_DEs_MathFuncs_Maple2015.15.2.mla.zip used by the whole system or only by the Physics package?

@Kitonum 

coeffs(a*x+b*y,[x,y]) =~ coeffs(d*x+e,[x,y]);

@Preben Alsholm 

Most probably yes. But solve cannot use your strategy I think.
Since dsolve gives the correct result for ODE with no IC, it is simpler to use:

res:=dsolve({ diff(y(t),t) = y(t)^2 - y(t)^3 });
eval(res,[t=0,y(t)=1/10]);
solve(%);
eval(res,%);

@Boby 

You must define an F if you want to see the derivatives. Otherwise F stands fo a generic function.

If you want F to be the polynomial   a x^2 + a x + b  , you must enter

F := x -> x^2 + a*x + b;

(in Maple it is called procedure).

After that you can simplify again your expresion (click in the for loop and press ENTER) and this polynomial will be derivated.

@Axel Vogt 

It is not constant as

plots[odeplot](solA0,[theta,r(theta)],theta=varepsilon..5/1000);

shows. It is only a precision problem, abserr and relerr are used probably with some heuristic tests which are not enough when the solution is singular.

If you hope for some help, please state the complete problem, without asking to google for "research papers" with no references/links.

@tomleslie 

Is seems to be hard to obtain o prescribed precision.

In your worksheet one obtains the value of A:

A0 := 199.961318896644828500644760320

and the solution

Now, if we increase a bit the precision in dsolve without parameters for this A0:

solA0 := dsolve( {Ode1, r(varepsilon) = A0, (D(r))(varepsilon) = 0},
               numeric,
               range = varepsilon .. 5*(1/1000),
               abserr=1.0e-14,
               relerr=1.0e-13
             );

one obtains for solA0(5/1000);

i.e.  r(theta) has only 8 correct digits!

How many digits can we guarantee for A0?

@tomleslie

I do not think that the optimization is essential here; fsolve will do it after adjusting Digits.

The real problem is the accuracy of dsolve with parameters.

@Preben Alsholm 

I think we should have serious concerns about the precision because of the singularity.

For example fsolve (instead of NLPSolve) refuses to solve the equation even in the interval  199.96..199.97.

If we set Digits:=27 before fsolve ==>

199.961306129991433460330154

but the problem is how good is the approximation given by dsolve with parameters.

@acer 

I did not do it on purpose.

showstat(plottools:-`transform/object`);

simply does not include it.

Anyway I don't think that maplesoft will have something to lose, on the contrary!