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]

I do not understand the other half of your code and wrote the animation in 3D anew:

with(plots):
Repltlist := [[sin(t), 0], [cos(t), (1/3)*Pi], [cos(2*t), (2/3)*Pi]]:
Impltlist := [[-sin(3*t), 0], [-2*cos(t), (1/3)*Pi], [-cos(4*t), (2/3)*Pi]]:
L1:=map(s->[s[1]*cos(s[2]), s[1]*sin(s[2]), t], Repltlist);
L2:=map(s->[s[1]*cos(s[2]), s[1]*sin(s[2]), t], Impltlist);
A:=a->spacecurve(L1, t=-Pi..a, color=red, thickness=3):
B:=a->spacecurve(L2, t=-Pi..a, color=blue, thickness=3):
animate(display, ['A'(a), 'B'(a), axes=normal, labels=[x,y,t] ], a=-Pi..Pi, frames=90);

The function  delta  is incorrectly defined. If  k=a  then  k-a=0  and  delta can be defined as follows:

delta := x->`if`(x=0, 1, 0);

However, Y[2], ... , Y[12]  can not be found by your for loop because it already for the first value  k=0  we have

k:=0:
add(delta(i-1)*(k-i+1)*(k-i+2)*Y[k-i+2], i = 0 .. k)+add((delta(i)-delta(i-1))*(k-i+1)*Y[k-i+1], i = 0 .. k)+lambda*Y[k] = 0;

                                                              A*lambda+B = 0

and this equation does not contain  Y[2] .
 

Edit. May be should be  delta(i)  instead of   delta(i-1)  in the first  add . Then everything works:

restart;
Y[0] := A; Y[1] := B;
delta:=x->`if`(x=0, 1, 0);
for k from 0 to 10 do Y[k+2] := solve(add(delta(i)*(k-i+1)*(k-i+2)*Y[k-i+2], i = 0 .. k)+add((delta(i)-delta(i-1))*(k-i+1)*Y[k-i+1], i = 0 .. k)+lambda*Y[k] = 0, Y[k+2]) end do;
y := sum(Y[j]*x^j, j = 0 .. 10);
 

At first I replaced the names  eq1[k_] .. eq4[k_]  with the simplier names  eq1 .. eq4  and I also made some more simplifications of your code.  The result is a system of 8 equations with 9 unknowns.  The  fsolve  command only solves a system in which the number of equations equals the number of unknowns. The solution can be obtained if one of the unknowns to give some value:

Download fdm-maple1.mw

If you want to move a circular vortex itself from right to left, the original system should depend on the parameter animation:

restart;
with(DEtools):
rho := 0.1:
w0 := 2:
sys := a->[diff(x(t),t) = y(t),diff(y(t),t) = -2*rho*y(t) - w0^2*(x(t)+a)];
P:=a->DEplot(sys(a), [x(t),y(t)], t = 0 .. 20-2*a, x=-2..2, y=-1.9..1.7, [[x(0) = cos(a)-a, y(0) = sin(a)]], scene = [x(t),y(t)], linecolor=blue, numpoints=1000):
plots:-animate(plots:-display,['P'(a), size=[600,300]], a=-0.7..1.4, frames=90);
                      

Edit.

The animation shows how the trajectory changes as the point  [x(0), y(0)] (the initial condition) moves along a semicircle from right to left:

restart;
with(DEtools):
rho := 0.1:
w0 := 2:
sys := [diff(x(t),t) = y(t),diff(y(t),t) = -2*rho*y(t) - w0^2*x(t)];
A:=plot([cos(s), sin(s), s=0..Pi], thickness=0):
P:=a->DEplot(sys, [x(t),y(t)], t = 0 .. 20,x=-2..2, y=-2..2, [[x(0) = cos(a), y(0) = sin(a)]], scene = [x(t),y(t)], linecolor=blue, numpoints=1000):
plots:-animate(plots:-display,['P'(a)], a=0..Pi, background=A, frames=90, scaling=constrained);

This is easy to do, only tickmarks on the y-axis you'll have to code yourself.

Example:

P:=plot([sin(x), cos(x)+2.4, 1.2], x=0..2*Pi, -1.2..3.6, color=[red,blue,black], tickmarks = [piticks, [seq(0.5*i = 0.5*i, i = -2 .. 2), seq(0.5*i+2.4 = 0.5*i, i = -2 .. 2)]], axes=box):
T:=plots:-textplot([[4.2,0.2,y=sin(x)], [4,2.5,y=cos(x)]], font=[times,roman,14]):
plots:-display(P, T, scaling=constrained);

Edit.

Adjustment:

restart;
Profit := (p-c)*d;
d:=a*p^(-b);
print(Solution); 
dProfit1 := diff(Profit, p); 
dProfit2 := diff(dProfit1, p); 
ddProfit2 := diff(Profit, p, p); 
pOpt := solve(diff(Profit, p), p);

Edit.