@nk2016

restart;
A:=-(cos((t/n)+t))/(sin((t/n)+t)):
sol:=solve(denom(A)=0, t, allsolutions):
indets(%)[1]:
h:=subs(%=1, sol):
r:=cos(t/n)^n:  r2:=eval(r, n=4):
x:=r2*cos(t):  y:=r2*sin(t):  
H:=seq(eval([x,y],n=4), t=[seq(h*i,i=0..4)]):
PP:=seq(plots:-display(plottools:-disk(H[i],`if`(i in [1,2,5],0.015,0.001), color=red)), i=1..5):
plots:-display(<plots:-display(PP, plot(r2, t=0..4*Pi, coords=polar, color=blue), title="The points at which the tangent is parallel to the axis Oy")| plots:-display(PP[3],PP[4], plot(r2, t=0..4*Pi, coords=polar, color=blue), view=[-0.05..0.05,-0.05..0.05], title="Zoom 20:1 near the origin")>, scaling=constrained);

Curves.mw


 

@nk2016   The addition to the previous code:

H:=seq(eval([x,y],n=4),t=[seq(h*i,i=0..4)]):
PP:=seq(plots:-display(plottools:-disk(H[i],`if`(i in [1,2,5],0.015,0.001),color=red)), i=1..5):
plots:-display(<plots:-display(PP, plot(r2, t=0..4*Pi, coords=polar, color=blue), title="The points at which the tangent is parallel to the axis Oy")| plots:-display(PP[3],PP[4], plot(r2, t=0..4*Pi, coords=polar, color=blue), view=[-0.05..0.05,-0.05..0.05], title="Zoom 20:1 near the origin")>, scaling=constrained);            

@asa12   

Prefix1:=proc(Expr)
if type(Expr, `*`) then return func1(Prefix1(op(1,Expr)), Prefix1(`*`(op(2..-1,Expr)))) else
if type(Expr, `+`) then return func2(Prefix1(op(1,Expr)), Prefix1(`+`(op(2..-1,Expr)))) else
if type(Expr, `^`) then return func3(Prefix1~([op(Expr)])[ ]) else
Expr  fi; fi; fi;
end proc:

Example of use:

Prefix1(a+b+c+d+e*c*d);
                               func2(func1(e,func1(c,d)),func2(a,func2(b,func2(c,d))))

@sumitzanje  It is easy to prove, that if we consider a[0], a[1], a[2], a[3], omega[s], L  as arbitrary constants, then your system is inconsistent. 

In fact, if we solve the system without the first equation, we get 2 solutions. If we now substitute these solutions in the first equation, then we do not get any identity:

restart;
Sys:={a[1]=sqrt((4*mu)/(3*alpha*omega^2)),a[2]=(2*beta*mu)/(9*alpha*omega^2),a[3]=sqrt(mu^3/(432*alpha*omega^2)), omega[s]=omega-(mu^2/16*omega)-((2*beta^2*mu)/(9*alpha*omega))};
Sol:=solve(Sys, {alpha,beta,mu,omega});

eval(a[0]=L+(2*beta*mu)/(3*alpha*omega^2), Sol[1]);
eval(a[0]=L+(2*beta*mu)/(3*alpha*omega^2), Sol[2]);

@dharr  I added 2 more ways to my answer and also uploaded the worksheet.

@nk2016   You do not have to type anything. Just copy it as text into the window of your worksheet and press the "Enter" key.

@Zeineb   if  F=1

@spalinowy  Why copy manually? In the file the simplified system has the name  Sys, and it can be used for further work with the system using the desired packages.

@spalinowy   I removed all called packages, including MTM, made restart at the beginning of the worksheet. Now it works.

Model_Maple3.mw

@spalinowy  Put here your system, or the link to the worksheet.

Should be  

MTM:-solve(5*x-3 = 19);

@vv   Very witty and beautiful solution, vote up! 

It's a pity that the symmetry is broken for any even N.

@vv  Problem point is that each line should intersect with each line. Only in this way we get the maximum number of pieces. I do not know whether there is such a turn-symmetric structure if  n>5

@Mac Dude  Of course this is not the completed code, but only the beginning. OP wants to solve the original partial differential equation by difference method and, therefore, replaces each derivative by the corresponding difference quotient.

@jga  It may be helpful for you that your function can be decomposed into partial fractions:

f:=1/(1+q)^4/(1-q^2)^2*q^4*(q^2+1)*(((q^6+q^5-q^4-2*q^3+q+1)^6/q^13+(q^6+q^5-2*q^3-q^2+q+1)^6/q^25)):
convert(f, parfrac);