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.

I made some corrections to your document, and now everything works.

system_solve1.mw

The  createDegree  procedure solves your problem for a general case:

restart;
createDegree:=proc(L::{set,list}, n)
local L1, Degrees, k, Degs, NestedSeq;
L1:=convert(select(t->degree(t)<=n, L), list);
Degrees:=degree~(L1);
k:=nops(Degrees);
Degs:=[seq(floor(n/Degrees[i]), i=1..k)];
NestedSeq:=proc(Expr::uneval, L::list)
local S;
eval(subs(S=seq, foldl(S, Expr, op(L))));
end proc;
[NestedSeq(`if`(`+`(seq(Degrees[j]*(i||j), j=1..k))<=n,`*`(seq(L1[j]^(i||j), j=1..k)), NULL), [seq(i||j=0..Degs[j], j=1..k)])][2..-1];
end proc:

Example of use:
L := {a-2*b, b^2, (a+c)^2}:
for n from 1 to 4 do
createDegree(L,n);
od;

This can be done in several ways. The easiest way is

restart;
plots:-pointplot3d([[1,3,5],[10,30,55],[50,70,25]], style=line, color=red, thickness=2);
 

Example:

sols := [[0,0], [6,0], [6,4], [0,4], [3,6], [6,4], [0,0], [0,4], [6,0]]:
plots:-display(seq(plot(sols[1..n], color=red,thickness=3)$5, n=1..nops(sols)), insequence=true);

Int(Pi*(2*ln(y+sqrt(y^2-4))-2*ln(2))^2, y=2..exp(1)+exp(-1)):
IntegrationTools[Change](%, y+sqrt(y^2-4)=t);  # prompt 1
simplify(value(%));
eval(%, exp(4)-2*exp(2)+1=(exp(2)-1)^2);   # prompt 2
simplify(%) assuming exp(2)-1>0;   # The final result
                         

                                           4*Pi*(exp(1)-4+5*exp(-1))

Addition. Unfortunately Maple itself could not without prompting simplify this:

simplify(sqrt(exp(4)-2*exp(2)+1));   # Should be  exp(2) - 1

Edit.

Here is a simple example with the same error:

restart;
f:=x->diff(g(x), x):
f(0);

     Error, (in f) invalid input: diff received 0, which is not valid for its 2nd argument
 

Maple just substitutes  x=0  and we have  diff(g(0), 0)

A workaround  to avoid the error is usage  D - operator instead of  diff :

restart;
g :=  (x,y)->D[1,1](u1)(x, y)+D[1,2](u2)(x, y);
f :=  y ->g(0, y);
f(0); 
f(1.2);

I understood that the word "simultaneously" means that in every moment  a = b = c :

restart;
eq1:=((diff(f(x),x$3)))+f(x)*diff(f(x),x$2)-a*diff(f(x),x$1)^2=0:
eq2:=(diff(g(x),x$2))+b*f(x)*diff(g(x),x$1)=0:
bc1:=f(0)=0,D(f)(0)=1,D(f)(5)=0,g(0)=0.5,g(5)=0:
sol:=dsolve(subs(a=0.5,b=0.5,{bc1,eq1,eq2}), numeric):
L:=evalf[5]([seq(eval(a*f(x)+b*diff(g(x),x)+c*diff(f(x),x)*g(x),[op(sol(1)),op([a,b,c]=~0.2*i)]), i=0..5)]):
Matrix([[a,b,c,Expr],seq([(0.2*i) $3,L[i+1]], i=0..5)])^%T;

Unfortunately an accurate visualization requires quite a hard work. Here is an example of the visualization of the body of the rotation  a parabola  around the x-axis by disc (washer) method:

f:=x->sqrt(4-x):
with(plots): with(plottools):
A:=plot(f, 0..4, color=black, thickness=2):
B:=plot(f, 0..4, color=yellow, filled=true):
C:=rectangle([2.4,f(2.5)], [2.6,0],color=blue):  # Area element
T:=textplot([[2.5,0,"x",align=below], [2.8,0.7,"f(x)"], [2.5,1.4,dx]], font=[times,roman,14]):
display(C,B,A,T, scaling=constrained, view=[-0.7..4.7,-0.7..2.7]);
A1:=plot3d([f(x)*cos(phi),x,f(x)*sin(phi)], phi=Pi/2..2*Pi, x=0..4, style=surface, color=yellow):
A2:=plot3d([r*cos(phi),0,r*sin(phi)], phi=Pi/2..2*Pi, r=0..2, style=surface, color=yellow):
A3:=plot3d([x,y,0], y=0..4, x=0..f(y), style=surface, color=yellow):
A4:=plot3d([0,y,z], y=0..4, z=0..f(y), style=surface, color=yellow):
B1:=plot3d([f(2.5)*cos(phi),x,f(2.5)*sin(phi)], phi=0..2*Pi, x=2.4..2.6, style=surface, color=blue):
C1:=plot3d([r*cos(phi),2.4,r*sin(phi)], r=0..f(2.5), phi=0..2*Pi, style=surface, color=blue):
C2:=plot3d([r*cos(phi),2.6,r*sin(phi)], r=0..f(2.5), phi=0..2*Pi,style=surface, color=blue):
display(A1,A2,A3,A4,B1,C1,C2, axes=normal, scaling=constrained);  # The body with a cutout for better visibility
Int(Pi*``(f(x))^2, x=0..4)=int(Pi*f(x)^2, x=0..4);  # Calculation of the volume