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

We have one equation with three unknowns in the real domain. There is an obvious trivial solution (0, 0, 0). Since the equation is homogeneous, if there is a nontrivial solution  (a, b, c) , then  (t*a, t*b, t*c)  also be a solution for any  t>0 . Therefore, the set of all solutions is a cone in R^3 . It is enough to examine the structure of the solutions on the unit sphere. Below we can see from the plot that at the unit sphere we have 3  isolated branches solutions. One solution has been found numerically.

restart;
Eq:=a^2+b^2+c^2+a*b+a*c+b*c-  (a+b-c)*sqrt(2*a*b+a*c+b*c)-(a+c-b)*sqrt(a*b+2*a*c+b*c)-(b+c-a)*sqrt(a*b+a*c+2*b*c):
Eq1 := eval(Eq, [a = cos(phi)*sin(theta), b = sin(phi)*sin(theta), c = cos(theta)]):
plots[implicitplot](Eq1, phi = 0 .. 2*Pi, theta = 0 .. Pi, color = red, grid = [1000, 1000]);
RootFinding[Analytic](eval(Eq1, theta = 0.8), re = 1..2, im = -0.1..0.1);  # a root
evalf(eval(Eq1, [phi = %, theta = 0.8]));  # Check