DuplicationMatrix:=proc(n)
local A, a, L1, L2;
uses ListTools;
A:=Matrix(n, symbol=a, shape=symmetric);
L1:=[seq(seq(A[i,j], j=1..n), i=1..n)];
L2:=convert(L1, set);
Matrix(n^2,n*(n+1)/2, {seq(seq((i,j)=`if`(j=Search(L1[i], L2),1, 0), j=1..n*(n+1)/2), i=1..n^2)});
end proc:

Example of use:
C:=DuplicationMatrix(3);

Check.

A, V1, V2:=<a,b,c; b,d,e; c,e,f>, <a,b,c,d,e,f>, <seq(seq(A[i,j], j=1..3), i=1..3)>;
C.V1=V2;
 

Under the condition that  u, v, delta, lambda, phi, beta, alpha  are in the interval  ]0,1]  there are no solutions.

The proof.

Let   A=(alpha+lambda+u+delta)/(u+v+alpha) ,  B=(alpha+lambda+u+delta)/(u+v)

From the first inequality we have  beta/phi>B  and from the fourth inequality we have  beta/phi<A .  We have a contradiction because  A<B .

Theta:=[ 0, Pi/3, 2*Pi/3, Pi]:
Phi:=[ 0, Pi/3, 2*Pi/3, Pi, 4*Pi/3, 5*Pi/3, 2*Pi, 7*Pi/3]:
F:=(theta,phi)->[cos(phi)*sin(theta), sin(phi)*sin(theta), cos(theta)]:
{seq(seq(F(theta,phi), phi=Phi), theta=Theta)};
A:=plots[pointplot3d](%, color=red, symbol=solidcircle, symbolsize=20):
B:=plottools[sphere]([0,0,0], 1):
plots[display](A, B);

a := [1, 4, 2, 6, 8]:  b:={a[ ]};
{seq(seq(parse(cat(i,j)), j=b minus {i}), i=b)};
select(`<`, b union %, 50);

                 {1, 2, 4, 6, 8, 12, 14, 16, 18, 21, 24, 26, 28, 41, 42, 46, 48}

or using  combinat:-permute  instead of nested  seq  :

a := [1, 4, 2, 6, 8]:  b:={a[ ]};
map(parse@cat@op, combinat:-permute(b,2));
select(`<`, b union {%[ ]}, 50);

Edit.
                     

Apply  evalf  command to your answer.

restart;
with(Student[NumericalAnalysis]):
f:=x->abs(x):
a, b:=-1., 1.:
h:=(b-a)/4:
xx:=[a,a+h,a+2*h,a+3*h,b]:
yy:=map(f, xx):
xy:=zip((u,v)->[u,v], xx, yy);
s2:=CubicSpline(xy, independentvar=x);
P:=expand(Interpolant(s2));
plot([f(x), P], x=-1..1);
Draw(s2);
M:=maximize(abs(abs(x) - P), x=-1..1);

Edit.

restart; 
1 = int(abs(A*exp(-(1/3)*x^2))^2, x = -infinity .. infinity); 
solve({%, A > 0}); 
A := eval(``(A^4)^(1/4), %);   # The final result

Addition. If in the future you are going to do any numerical manipulations with  A , you must first use  expand  command, for example:

evalf(expand(A));

                                  0.6787185469

Examples:

L:=[x1+x2-x3=0, 2*x1-x2+x3=0, x2-3*x3=0]:
convert(L, Vector);
# or
%piecewise(seq(op([``, L[i] ]), i=1..nops(L)));

Addition. You can also display it in a column without any brackets using a for loop:

L:=[x1+x2-x3=0, 2*x1-x2+x3=0, x2-3*x3=0]:
for i from 1 to nops(L) do
L[i];
od;
 

I think the reason is incorrect conversion, which does not need:

int(f(x)*f(y)*x^2*abs(x+y), x=-infinity..infinity, y=-infinity..infinity);
evalf(%);

                                     3/sqrt(Pi)

                                 1.692568750

f:=(x,n)->(1-x^2)^n/n!:
is(diff(f(x,n),x,x) = 2*(1-2*n)*f(x,n-1)+4*f(x,n-2));

                                     true
 

The proof below uses only the initial Maple conversion and  simplest transformations with roots:

Download 111.mw

I  simplified a little notations for your system. Maple finds all solutions of the system in a very bulky form. To simplify the answer the numerical values of the system parameters should be specified. See the attached file.

solutions1.mw

Edit.

Transformation to the  simple_answer :

ode:=diff(y(x),x) = sec(x)^2*sec(y(x))^3:
dsolve(ode, implicit);
A:=3*subs(sin(x)/cos(x)=tan(x), expand(-%));  # An intermediate result
factor(subs(_C1=_C1/3, A-add(op(i, lhs(A)), i=[1,4])));   # The final result

Edit.

Try to replace the code by this:

dsys := {2*m1*(a+l*sin(phi(t)))^2*(diff(diff(theta(t), t), t))+4*m1*(a+l*sin(phi(t)))*l*cos(phi(t))*(diff(theta(t), t))*(diff(phi(t), t)) = M, 2*m1*l^2*(diff(diff(phi(t), t), t))+4*m2*l^2*sin(2*phi(t))*(diff(phi(t), t))*(diff(phi(t), t))+4*m2*l^2*sin(phi(t))^2*(diff(diff(phi(t), t), t))-2*m1*(a+l*sin(phi(t)))*l*cos(phi(t))*(diff(theta(t), t))*(diff(theta(t), t))-2*m2*l^2*(sin(2*phi(t)))*(diff(phi(t), t))*(diff(phi(t), t)) = -(2*(m1+m2))*g*l*sin(phi(t))-2*k*l^2*sin(2*phi(t)), phi(0) = 0, theta(0) = 0, (D(phi))(0) = 0, (D(theta))(0) = 0};

{M = 10, a = .5, g = 9.81, k = .1, l = .5, m1 = 10, m2 = 1}:
dsn1 := dsolve(eval(dsys, %), numeric);

Edit.  Multiplication sign was missed. See

You can easily automate your manual solution using  minimize  command:

restart;
y0:=solve(x^2*y = 16, y);
minimize(eval(4*x^2+4*x*y, y=y0), x=0..infinity, location);
x0:=op([1,1,1], %[2]);
Minimum=%%[1], [op([1,1,1],%%[2]),y=eval(y0,x0)];  # The final result

                              Minimum = 48, [x = 2, y = 4]