Simplest way, I think, is

> fprintf("D:/Data/Data6/My_File7.txt", "%f", Q);

 as suggested by the error message.

First, I would avoid the RealDomain package, though I don't think it makes a difference in this case.  Here's the meat of the issue.

> Q:= convert(P_HAN_s_dif,rational);

Q := (1/2)*(3000*Pv+10000*T+5002*PLu-2500*Pu)/sqrt((30*Pv+100*T+50*PLu-25*Pu)^2+PLu^2)+(1/2)*(10000*T+5002*PLu+2900*Pv-2502*Pu) /sqrt((-100*T-50*PLu-29*Pv+25*Pu)^2+(Pu-PLu)^2)

> S := [solve(Q, PLu)];

S := [-(5/2501)*(200000*T^2-50020*T*Pu+118000*T*Pv+5*Pu^2-14506*Pv*Pu+17400*Pv^2)/(200*T+59*Pv), 5*Pu*(6*Pv-5*Pu+20*T)/(Pv-50*Pu)]

Now the usual way of solving equations involving radicals is to convert them to polynomials and solve the polynomials.  Unfortunately, the solutions you get are not always solutions of the original equations, because sometimes they may require taking the wrong branch of the radicals.   If you want to filter out these spurious solutions, you should substitute the solutions into the original equation and check if they work.  Unfortunately, Maple doesn't always do this.  In our case:

> V1:=simplify(eval(Q,PLu=S[1])) assuming real;

V1 := -(-signum(29*Pv+100*T+25*Pu)*Pu+Pu*signum(-6*Pv+5*Pu-20*T)-50*Pv*signum(29*Pv+100*T+25*Pu)+50*signum(-6*Pv+5*Pu-20*T)*Pv)* signum(29*Pv+100*T+25*Pu)*signum(-Pv+50*Pu)*signum(-6*Pv+5*Pu-20*T)/sqrt(Pv^2+Pu^2)

> V2 := simplify(eval(Q,PLu=S[2])) assuming real;

V2 := -sqrt(2501)*(-50*Pv*signum(29*Pv+100*T+25*Pu)+50*Pv*signum(6*Pv-5*Pu+20*T)-signum(29*Pv+100*T+25*Pu)*Pu+ signum(6*Pv-5*Pu+20*T)*Pu)*signum(29*Pv+100*T+25*Pu)*signum(200*T+59*Pv)*signum(6*Pv-5*Pu+20*T) /sqrt(5981*Pv^2+23600*T*Pv+100*Pv*Pu+40000*T^2+Pu^2)

Now either of V1 and V2 could be 0 or nonzero, depending on the signs of various linear combinations of Pv, Pu and T.  Note that at this stage we haven't told Maple anything about these signs, so we can't blame Maple if one of them is nonzero for the particular values we have in mind.

 > eval([V1,V2], {T=-40.11096379,Pu=-50.84351320,Pv=-141.8108432});

[94.80773532, 0.]

If the vertex is at the origin, then the radius is not a function of theta and phi, because the surface consists of lines at the same theta and phi.  You could do it parametrically, e.g. for a cone of opening angle alpha (where alpha is a number from 0 to pi):

> plot3d([r, theta, alpha], r=0..1, theta = 0 .. 2*Pi, 
     coords=spherical, scaling=constrained);

> de := diff(x(t),t) + (1-cos(t))*x(t)=exp(sin(t));

> dsolve(de);

x(t) = exp(sin(t))+exp(-t+sin(t))*_C1

See the help page ?dsolve for more options (there are many...)

ApproximateInt works correctly if you use an expression instead of a function, i.e.

ApproximateInt(f(t),t=x..y,method=simpson)

However, the syntax ApproximateInt(f, x..y) does not work.  In fact, the result does not depend on f at all.

> ApproximateInt(f, a .. b);

1/2*b^2-1/2*a^2

[ The <maple> tag is acting up again: that should be 1/2*b^2 - 1/2*a^2 ] 

What seems to be happening is that Maple interprets this input as

ApproximateInt(f, f=a..b)

To be fair, the help page ?ApproximateInt does not mention this syntax at all, so ApproximateInt(f, x..y) is an error on the part of the user.  However, this is definitely something that should be caught (either returning the same result as ApproximateInt(f(t), t=x..y) or an error message).  I'm submitting an SCR.

It's somewhat tricky, I think, because display_allGR displays things rather than returning a result.  But you can do things such as

> R[compts]:= convert(series(R[compts], omega, 3), polynom);

> with(LinearAlgebra):

> M:=Matrix(16, 16, [[1,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0] , [0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0], [1,0,1,0,0,0,0,0,0,0,0,0,1,0,1,0], [0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,0], [1,0,0,1,0,0,0,0,1,0,0,1,0,0,0,0], [0,0,0,0,1,0,1,0,0,0,0,0,1,0,1,0], [0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0], [0,0,0,0,1,0,0,1,1,0,0,1,0,0,0,0], [0,1,1,0,0,0,0,0,0,0,0,0,0,1,1,0], [1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1], [0,1,0,1,0,0,0,0,0,1,0,1,0,0,0,0], [0,1,0,1,0,0,0,0,0,0,0,0,0,1,0,1], [0,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0], [0,0,0,0,1,0,0,1,0,0,0,0,1,0,0,1], [0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0], [0,0,0,0,0,1,0,1,0,0,0,0,0,1,0,1]]);

> Rank(M);

              9

So M is a singular matrix with rank 9.  A basis of the null space of its transpose is:

> C:= NullSpace(M^%T);

So the system M.X = b is solvable if and only if c^T . b = 0 for all c in C.

Now your b is

> b:=<exp(2*B*J_1+B*J_2),1,1,1,1,exp(-2*B*J_1+B*J_2),
  exp(2*B*J_1-B*J_2),exp(-2*B*J_1-B*J_2),exp(-2*B*J_1-B*J_2),
 exp(2*B*J_1-B*J_2),exp(-2*B*J_1+B*J_2),1,1,1,1,exp(2*B*J_1+B*J_2)>;

> C[1]^%T . b;

4-2*exp(2*B*J_1-B*J_2)-2*exp(-2*B*J_1-B*J_2)

(and similarly for the others).  Your system won't have a solution unless that expression happens to be 0

This has nothing to do with being in a proc or not.  The code

(evalf@f)()/10^Digits 

returns a number (a different one each time it is evaluated).   In your first attempt,

seq((evalf@ f))()/10^Digits, i = 1 .. 5);

evaluates that five times, returning five different numbers.  In your second attempt,

U := uniform(1..3);

evaluates it once, returning one number which is assigned to U.  Now each time you call U(), you will get that number (a number can be considered as a function, namely the constant function whose value is the number).  If you tried instead

seq(uniform(1..3), i=1..5);

you would again get five different numbers.

This is a rather complicated nonlinear PDE, with no obvious way of separating variables.  Do you have any reason to think that it has  closed-form solutions?

Well, actually there are some: e.g.

> eq:= %;
   eq1:= eval(eq, f(theta, `&varphi;`) = g(theta) + K*`&varphi;`);
   dsolve(eq1);

g(theta) = Int(1/6*6^(1/2)/c/C*(cos(theta)^2*(-A^2*omega^2*cos(theta)^2+A^2*omega^2*cos(theta)^4+6*C^2*c^2*_C1-6*C^2*K^2*c^2-6*C^2*c^2*_C1*cos(theta)^2-6*K^2*ln(cos(theta)+1)*C^2*c^2+6*K^2*ln(-1+cos(theta))*C^2*c^2*cos(theta)^2+6*K^2*ln(cos(theta)+1)*C^2*c^2*cos(theta)^2-6*K^2*ln(-1+cos(theta))*C^2*c^2))^(1/2)/sin(theta)/cos(theta)^2,theta)+_C2

g(theta) = Int(-1/6*6^(1/2)/c/C*(cos(theta)^2*(-A^2*omega^2*cos(theta)^2+A^2*omega^2*cos(theta)^4+6*C^2*c^2*_C1-6*C^2*K^2*c^2-6*C^2*c^2*_C1*cos(theta)^2-6*K^2*ln(cos(theta)+1)*C^2*c^2+6*K^2*ln(-1+cos(theta))*C^2*c^2*cos(theta)^2+6*K^2*ln(cos(theta)+1)*C^2*c^2*cos(theta)^2-6*K^2*ln(-1+cos(theta))*C^2*c^2))^(1/2)/sin(theta)/cos(theta)^2,theta)+_C2

where the integral can be done in closed form in the case K=0.

The problem is that there is no (...)^(1/2) in the denominator, since your expression is actually represented as (numerator)*(...)^(-1/2).  One thing you can do is this.

> Q:= simplify(f(theta));

> S := indets(Q, sqrt);

S := {1/((a^2*cos(theta+1/4*Pi)^2+b^2-b^2*cos(theta+1/4*Pi)^2)^(1/2)), (a^2*cos(theta+1/4*Pi)^2+b^2-b^2*cos(theta+1/4*Pi)^2)^(1/2)}

[ By the way, that was produced by clicking the red maple-leaf icon in the editor: you copy the output from Maple and paste it into the "Maple Math Expression" field in the dialog that pops up]

> subs(S[1] = 1/tmp, S[2]=tmp, Q);

(a*b+c*cos(4*theta)*tmp)/tmp

With that initial condition, the animation is going to be rather boring, so I'll try a different one.

> de := diff(theta(t), `$`(t, 2))+c*(diff(theta(t), t))+9.8*sin(theta(t))/L = 0; c := 2; L := 2;
init := theta(0) = 0, (D(theta))(0) = 4;
sol := dsolve({de, init}, theta(t), numeric, output=listprocedure);
Theta := subs(sol, theta(t));
plots[animate](plots[arrow],[[0,0],[L*sin(Theta(t)),-L*cos(Theta(t))],colour=red],t=0..5, axes=none,scaling=constrained);

In a Maple plot, the idea is to put all the elements together at once rather than plotting "incrementally", one element at a time.  Thus your loop could produce sets of blue, green and red points, then you could use pointplot (in the plots package) to produce
plots of these points, and display (also in the plots package) to combine them all together.

That is a good question.  The approach that I would use is to take the differential equation with the first data point as initial condition, consider the values at the other t's in the data set as functions of the parameters, and use nonlinear least squares to fit these to the data.  The LSSolve command can use procedures (in particular, involving dsolve(..., numeric)).   It will need the derivatives of the x and y values with respect to the parameters: these can be obtained by enlarging our system of differential equations.  I'll write e.g. xa(t) as the derivative of x(t) (regarded as a function of the parameters as well as t) with respect to a.  Differential equations for these can be obtained by differentiating the original system. 

Note that there's an extra degree of freedom in the parameters: you could multiply all of a,b,c,d,f by some nonzero factor and divide m and n by it, without changing the equations.  So I'll assume m=1.
 

> sys := {D(y)(t) = a+b*y(t)+c*x(t)+d*x(t)*y(t)+f*y(t)^2,
    D(x)(t) = n*(a+b*x(t)+c*y(t)+d*x(t)*y(t)+f*x(t)^2)};
 

> param := [n,a,b,c,d,f];

Here are the equations for the derivatives with respect to parameters:

> Derivs:= {seq(seq(D(cat(v,p))(t) =  subs(
  diff(x(t,p),p)=cat(x,p)(t),diff(y(t,p),p)=cat(y,p)(t),
   x(t,p)=x(t), y(t,p)=y(t), diff(subs(x(t)=x(t,p),y(t)=y(t,p),  
  subs(sys,D(v)(t))),p)),p=param),v=[x,y])};

I'll try it with the following data:

> Data := [[0., 0., 1.], [1., .47, .94], [2., .73, .92], 
  [3., .86, .93], [4., .93, .95], [5., .96, .96]];
  N := nops(Data);

The initial conditions for x and y are given by the first data point; the derivatives
with respect to the parameters will be 0 there.

> ics := {x(0) = Data[1,2], y(0) = Data[1,3], 
     seq(seq(cat(v,p)(0) = 0, p = param), v = [x,y])};
 

Maple 12 provides a nicer way of dealing with parameters in dsolve(..., numeric).
Previous versions needed to use assignments to global variables.

> Sol := dsolve(sys union Derivs union ics, numeric, 
        parameters = param);

The procedure Res will provide the residuals, i.e. the differences between the values of x and y at the given t's and the values given in Data.  It is passed the parameter values as a Vector in the first argument, and returns the values of the residuals as a Vector in the second.

> Res:= proc(X,Y)
   local j, R;
   Sol(parameters = convert(X,list));
   for j from 2 to N do
      R:= Sol(Data[j,1]);
      Y[2*j-3]:= subs(R, x(t)) - Data[j,2];
      Y[2*j-2]:= subs(R,y(t)) - Data[j,3];
   end do
end proc;

The procedure Jac will provide the Jacobian matrix of the residuals with respect to the parameters.  It is passed the parameter values as a Vector in the first argument, and returns the values as a Matrix in the second.

 

> Jac := proc(X,Y)
    local j,k, R;
    Sol(parameters = convert(X,list));
    for j from 2 to N do
       R:= Sol(Data[j,1]);
       for k from 1 to 6 do
          Y[2*j-3, k] := subs(R, cat(x,param[k])(t));
          Y[2*j-2, k] := subs(R, cat(y,param[k])(t));
       end do
   end do
end proc;

And now for the optimization:

> with(Optimization):
   LSSolve([6,2*N-2],Res, objectivejacobian=Jac,
     initialpoint=<2,0.1,-0.1,0.2,-0.1,-0.1>);
 

[.10658455836617e-4, RTABLE(259413740,MATRIX([[1.47678922329157314], [-.909293770783568212], [1.10226087521962302], [1.31994271304865984], [-1.25914681017952136], [-.263174744443265829]]),Vector[column])]

To see the results for these parameter values:

> Sol(parameters = convert(%[2], list));
   Matrix([seq(eval([t, x(t), y(t)], Sol(Data[j,1])),j=1..N)]);

RTABLE(263739708,MATRIX([[0., 0., 1.], [1., .470531728907655433, .939867825078991136], [2., .728250451686391020, .920382998001185282], [3., .862609215916281924, .930397624753999453], [4., .928797484673319352, .947726261608177522], [5., .959503363012274788, .961995096046812659]]),Matrix)

> with(plots):
display(pointplot3d(Data,colour=black,symbol=box),
 odeplot(Sol,[t,x(t),y(t)],t=0..5),axes=box,     
 orientation=[-76,64]);

One problem with this approach is that for some values of the parameters we get error messages such as

Error, (in Optimization:-LSSolve) cannot evaluate the solution further right of 3.1217572, probably a singularity

presumably due to the fact that the solution "blows up".  This might be fixed by modifying the differential equations to keep the solutions bounded.

Your loop certainly does something: it used up all the memory on my computer, causing it to grind to a halt: I eventually had to use the Task Manager to shut it down after about 2256 iterations.  Something in your program seems to be eating up memory. 
 

Maple is quite correct in this case: x = 6 is the only solution (if, as Maple does, you use the principal branch of the square root).  The quadratic equation obtained from your equation does have two roots, x=6 and x=150, but x=150 is not a solution of your equation.