@John Fredsted  Thanks for this elegant and simple solution (thumb up)!

I propose an improvement. Now it can also split polynomials and radicals.

Splitting:=proc(Expr)
local x, y, Expr1, Expr2;
x, y:=indets(Expr, name)[];
Expr1, Expr2 := selectremove(z -> has(z,x) and has(z,y), simplify(factor(Expr), symbolic));
``(simplify(select(has, Expr2*expand(Expr1), x)))*
``(simplify(remove(has, Expr2*expand(Expr1), x)));
end proc:

Example of use:

Expr:=(x^2*y+y)*y^3*exp(-x-y+1)*3^(-x-y)*sqrt(x^2*y-y);
Splitting(Expr);
                    

@AndreaAlp  You wrote  "I'm talking about the original surfaces f(x,y) =z and g(x,z)=y. If I understood correctly, f is a plane and g is a parabolid. But I cannot find the correct centre of the intersecting ellipse."
In fact, these two surfaces do not intersect. Below we show it strictly analytically (although it is easy to see also from the plotting). First, we eliminate z-variable from the system, and then reduce the resulting equation to completed square form:

restart;
with(Student:-Precalculus):
my_plane := z=2024.30587449691-.341275505799078*x-3.89936179341114*y:
my_quadric := y=-10595.4104931095+6.73749956241827*x+42.1022818380012*z-0.654818649508000e-1*x^2-0.421174257681115e-1*z^2:
Sys:=[my_plane, my_quadric];
A:=eliminate(Sys,z);
B:=CompleteSquare(%[2][], x);
C:=op(1,B)+CompleteSquare(`+`(op(2..-1,B)), y)=0;  # The equation C has not any real solutions
plots:-implicitplot3d([my_plane, my_quadric], x=0..150, y=-100..450,z=400..600, style=surface, color=["LightBlue", yellow], numpoints=10000, scaling=constrained, axes=normal); # The visial check
                                

Let us dwell on the version with surfaces  z=f(x,y) and g(x,z)=0 . Of course, the first surface is a plane, and the second surface is an elliptical cylinder. To save all proportions, scaling=constrained  option is used:

restart:      
with(plots):    
with(Student:-Precalculus):
my_plane := z=2024.30587449691-.341275505799078*x-3.89936179341114*y:
my_quadric := -10595.4104931095+6.73749956241827*x+42.1022818380012*z-0.654818649508000e-1*x^2-0.421174257681115e-1*z^2:
Expr:=CompleteSquare(my_quadric, [x,z]);
x0:=solve(op(2,Expr))[1]; z0:=solve(op(1,Expr))[1];  # [x0,z0] are the coordinates the intersection of the axis of the cylinder with (x,z)-plane
A:=spacecurve([x0,t,z0], t=300..500, color=violet, thickness=3, linestyle=2):
B:=intersectplot(subs(y=y-1,my_plane), my_quadric, x=0..100, y=300..500,z=420..580, color=red, thickness=4): 
P:=plots:-implicitplot3d([my_plane, my_quadric], x=0..100, y=0..430,z=420..580, style=surface, color=["LightBlue", yellow], numpoints=10000, scaling=constrained, axes=normal, lightmodel=light4):
display(A,B,P);

Here is another solution, also based on reduction to the solution of a differential equation, but without any substitutions. The values of all the parameters for the numerical solution I took arbitrarily:

Download System.mw

Here is the corrected solution. I did not delete the previous version, because it is more simple and the idea of a solution is clearer on it:

System_new.mw

Edit.

You have a system of 5 linear ordinary differential equations with 14 parameters. Maple can not solve it symbolically, but it easily solves numerically. For this, we must specify the values of all parameters and the initial conditions. Here is an example of a solution with visualization: 
 

Download System.mw

restart;
RealDomain:-solve({x=r*cos(theta), y=r*sin(theta)}, {r,theta}, explicit);

Formally, both solutions are correct. But if you want  r>=0 , then leave only the first solution. Your expected result is incorrect, because  arctan(y/x)  returns the angle in the range  -Pi/2 .. Pi/2  only. See help on arctan  function about  two-argument function arctan(y, x)

restart; 
with(LinearAlgebra):
for n from 5 to 8 do:
for i from 1 to n do for j from 1 to n do A[i,j]:=0 od: od:
for i from 1 to n-1 do A[i,i+1]:=3 od:
for i from 1 to n do A[i,i]:=4 od:
for i from 2 to n do A[i,i-1]:=5 od:
A[n]:=Matrix([seq([seq(A[i,j], j=1..n)], i=1..n)]);
P[n]:=Determinant(evalm(1-t*A[n]));
od:
seq(P[n], n=5..8);

Explanation. Your mistake is that at each step of the loop on n-variable you use the same name  A , which also coincides with the name of the indexed variable  A .

Your code can be written much shorter:
restart;
with(LinearAlgebra):
A:=n->Matrix(n, (i,j)->piecewise(i=j,4, j=i+1,3, j=i-1,5, 0));
seq(A(n), n=5..8);
seq(Determinant(evalm(1-t*A(n))), n=5..8);
 

Edit.

restart;
with(plots):
m1:=70.0:  m2:=25.0:  m3:=3.0:  g:=9.8:
Sys:=[T+N1-m1*g=0, N2-N1-m2*g=0,T-N2-m3*g=0]:
solvec:=fsolve(Sys);
convert(evalf[4](solvec), list);  # Convert to the list and remove extra zeros
sol:=rhs~(%);  # Extract numeric values
A:=plots:-implicitplot3d(Sys, N1=0..300, N2=0..600, T=0..600, style=surface, color=[yellow, "LightBlue", green], axes=normal):  # Three planes
B:=pointplot3d(sol, color=red, symbol=solidsphere, symbolsize=18):  # The point in the space
display(A, B);   # Visualization of the system solution

Addition. If you want to get solutions in symbolic form as a function of parameters  m1 , m2 , m3  (as you wrote  "Actual true answers"), then you have at first to remove assignments from these parameters:

m1, m2, m3 := 'm1', 'm2', 'm3':
solvec:=solve({T+N1-m1=0, N2-N1-m2=0,T-N2-m3=0}, [T,N1,N2])[]; 

I think that the issue is in incorrectly setting the boundary conditions. Their number should be 8. Here everything is correct. But  D(f)(0)  and  D(g)(0)  are specified through  f(0)  and  g(0), which are unknown. If you directly specify  D(f)(0)  and  D(g)(0)  as  any numbers, everything will be OK, for example

bcs := h(0) = 0, p(0) = 0, theta(0) = 1, (D(f))(0) = 0.5, (D(g))(0) =  -0.5, f(etainf) = 0, g(etainf) = 0, theta(etainf) = 0; 

If we look closely at the marks on the t-axis, we see that all the graphs are constructed in a logarithmic scale. I did not understand why this happened, because I did not find the appropriate settings in your document. By default Maple uses a uniform scale on both axes. I just restarted your document and got other graphics that do not cause questions:

Download schechter_guo_v2_new.mw

In your code there are two "if", but only one "end if" .  Also what means  epsilon^() ?

The last initial condition  limit(diff((diff(f(x), x))/x, x), x = 0) = 0  follows from two conditions  D(f)(0)=0  and  (D@@2)(f)(0)=0 . But then the equation is obtained overdetermined (5 conditions instead of 4)  and Maple shows an error. The equation can be solved with the condition  D(f)(0)=0  and not very high accuracy:

restart;
ode := x^3*(diff(f(x), x, x, x, x))+alpha*(x^4*(diff(f(x), x, x))+x^3*(diff(f(x), x, x))-x^2*(diff(f(x), x)))-2*x^2*(diff(f(x), x, x, x))+3*x*(diff(f(x), x, x))-3*(diff(f(x), x))+R*x*f(x)^2-R*x^2*(diff(f(x), x, x))*(diff(f(x), x))-3*R*x*f(x)*(diff(f(x), x, x))+3*R*f(x)*(diff(f(x), x))+x^2*R*f(x)*(diff(f(x), x, x, x))-M^2*(x^3*(diff(f(x), x, x, x))-x^2*(diff(f(x), x))) = 0:
ics := f(0) = 0, f(1) = 1, D(f)(1) = 0, D(f)(0)=0:
sol:=dsolve(eval({ode, ics},[alpha=1,R=2,M=3]), method=bvp[midrich],numeric, maxmesh=1000, abserr=10^(-3));
plots:-odeplot(sol, [x,f(x)], x=0..1, color=red); 

Edit.

f:=(R,h,l)->[(1/3)*Pi*R^2*h, Pi*R*l, Pi*R^2+Pi*R*l]:

Examples of use:

f(3,4,5);
map(t->f(t[]), L);
                          

You do not need any loops. See this toy example:

A:=Matrix([[["BrakePower", "BrakePower"], [2.356, 2.356], [2.749, 2.749], [3.142, 3.142], [3.534, 3.534], [3.927, 3.927], [4.32, 4.32]], [["BrakePower", "S2"], [2.356, .303], [2.749, .271], [3.142, .256], [3.534, .249], [3.927, .244], [4.32, .241]], [["BrakePower", "S4"], [2.356, .256], [2.749, .225], [3.142, .211], [3.534, .205], [3.927, .2], [4.32, .197]]]):
y := a*x^2+b*x+c:
CurveFitting[LeastSquares](convert(A[1,2..-1], list), x, curve = y);
 

For your "polynomial" you can do so:

norm(subs(1.25=125, aa), infinity);

Here is a more general method that is suitable for a "polynomial" with any number of terms with fractional exponents:

max(abs~(map(p->`if`(type(evalf(op(1,p)),numeric),op(1,p),1),[op(aa)]))[]);
 

Edit.

Should be  s:=s+i  instead of  s:=s+a[i]

Doubt_on_list_new.mw

Explanation:  i  is a member of the list  a , and  a[i]  is the member of the list  a standing at  i_th place.