But you calculate completely different integrals, so the results are different.

V = Int(N, T__2 = 0 .. infinity)<>IntegrationTools:-Parts(V, exp(sigma-delta*(T__2-T__F)-(T__2-T__F)^2/T__p^2), T__2)

Using the option  T__2  you change the source integral. See the help on the  IntegrationTools:-Parts  command.

Take a look at a simple example in which we see that adding a third parameter in the  IntegrationTools:-Parts  command completely changes the original integral:

V:=Int(x^2*sin(x), x);
IntegrationTools:-Parts(V,x^2);
IntegrationTools:-Parts(V,x^2,x);

Your surface is the part of the green surface above the blue plane:

f:=x^3*y*z;
plot3d([solve(f=0.5, z), y], x=0..y, y=0..3, color=[green,blue], view=0..3, axes=normal);

The lowest point of your surface has equal coordinates  0.5^(1/5)

As alternative you can also use the  simplify  command with side relations option:

eqn:= N*x*y - x*y + f(x,y) = 0;
simplify(eqn, {x*y= 1});

This can be done in many ways. Here is a solution using a procedure that I wrote many years ago. All comments in the procedure are written in Russian. The full text of the procedure itself can be found in the attached file, and here its application to your example:

Download QuadricCurveAnalysis.mw

This is a known bug in  RealDomain  package when it does not take into account the domain of the equation. Here is a simple workaround:

restart;
RealDomain:-solve({(x-1)*sqrt(x^2-4) = 0, x^2-4>=0}, x);

                                     {x = -2}, {x = 2}

restart;
f:=t->piecewise(t<1 and t>=-1,t/2+1/2,t>=1 and t<2,-t+2);
F:=t->f(t-floor((t+1)/3)*3); #  F  is a periodic extension of  f  to the whole real axis
F(4.5);
plot(F(t), t=-4..5, scaling=constrained, size=[1200,150]);

Use the  LinearAlgebra:-GenerateMatrix  command for this (3 steps):
 

Download MatrixForm.mw

If you do not like the order of the elements in the set, then you can convert it into the list, and then sort this list according to the desired order:

E:={{1,2,3,4,5,6},{4,5,6,7,8,9,10,11},{10,11,12,13,14}};
        E := {{10, 11, 12, 13, 14}, {1, 2, 3, 4, 5, 6}, {4, 5, 6, 7, 8, 9, 10, 11}}
E1:=convert(E, list);
       E1 := [{10, 11, 12, 13, 14}, {1, 2, 3, 4, 5, 6}, {4, 5, 6, 7, 8, 9, 10, 11}]
sort(E1, (x,y)->x[1]<y[1]);
       [{1, 2, 3, 4, 5, 6}, {4, 5, 6, 7, 8, 9, 10, 11}, {10, 11, 12, 13, 14}]

Of course, to avoid this hassle, you should immediately specify E as a list and not as a set.

See help on the  diff  command.

The code for your example:

w:=exp(-beta*t-mu*x)*u(x,t);
diff(w, t);

Eq:=v1(t) - R1*(-i3(t)/n13 - i2(t)/n12) - L1*diff(-i3(t)/n13 - i2(t)/n12, t) = n12*(v2(t) - R2*i2(t) - L2*diff(i2(t), t));
Eq1:=expand((lhs-rhs)(Eq));
Terms:=[op(Eq1)];
DTerms:=select(has, Terms, diff);
Lhs:=add(DTerms);
collect(Lhs, diff)=collect(-(Eq1-Lhs), R1);  # The final result

You must specify the integration range, but you cannot use the symbol to as a name because it is a reserved word. I replaced it with  a :

a:=3:
evalf(int(exp((5+I*6)*x)*sin(1+erf(x)),x=-a..a));
evalf(int(exp((5+I*6)*x)*cos(1+erf(x)),x=-a..a));

restart;
eqn1 := e1*i1 + e2*i2 + e3*i3 = 0;
A:=i1=expand(solve(eqn1, i1));
subs(e2=e1/n12, A);

Or more automatically:

restart;
eqn1 := e1*i1 + e2*i2 + e3*i3 = 0;
A:=i1=expand(solve(eqn1, i1));
eqn2 := n12 = e1/e2;
subs(e2=solve(eqn2, e2), A);

Or using the  applyop  command:

restart;
eqn1 := e1*i1 + e2*i2 + e3*i3 = 0;
A:=i1=expand(solve(eqn1, i1));
applyop(simplify, [2,1], A, {e1/e2=n12});

Edit.
 

Use  plots:-pointplot3d  and  plots:-display  commands for this.

An example:

pointplot3d.mw

You forgot to specify values for  xbas, tbas :

Download anim_plot3d.mw
 

I rewrote your code as 2 procedures. The first one generates a full random polynomial of degree  N  with coefficients from the range  R . The second finds all rational roots of a polynomial with integer coefficients:

Download rational_zeros_new.mw