@one man After the transformation x^3 = s, the roots within each cluster are pretty much equally spaced, so it would be difficult to miss any one of them.

There is, however, the following consideration.  Suppose that we are looking for the roots of the function f(x) = x^2 + e.  We see that there are no roots if e > 0 and two roots if e < 0.  But if the size of e is smaller than what the floating point representation can resolve, deciding whether e is positive or negative is ill-posed, and that can lead to lost roots or false roots.

That same thing can happen in the problem that you have posed.  To detect that, I we would do the calculation with some setting of Digits, and then repeat with a larger value of Digits and check if the number of roots changes.

But you are asking for an a priori guarantee for detecting such exceptional cases.  I don't see a good way of doing that.

@one man The number of roots within each cluster is proportional to the length of the interval between the consecutive zeros of the function cos(0.25*surd(s,3)).  These intervals become longer as s grows, and so the cluster sizes become larger.  The first cluster, located in 1800 < s < 2181, is the smallest one, and as we have seen, it contains 122 roots.  In the same way we find that the next cluster corresponds to 15128 < s < 16646 and contains 482 roots.

As the clusters grow in size the farther away we go from the origin, I don't understand what you mean by the "most distant" cluster.

@nm Yes, I agree.   If Maple cannot determine the arbitrary constant, it should return NULL instead of something involving the constant.

@acer And it turns out that the "simplify" is not necessary:

@KIRAN SAJJAN Maple's pdsolve() is not the right tool for solving that problem.  I wrote earlier that pdsolve() is designed for solving evolution equations (initial value problems) while what you need is software for solving a boundary value problem. Maple's numerical PDE solver is not designed for that.

Additionally, I noted earlier that you need to take the outflow condition—equation (3) in the paper—into account.  That's missing in your worksheet.  That condition is not optional; the problem cannot be solved without it.

I am afraid that there is no easy solution to this problem.  If you really need to solve it, you will need to write your own finite elements (FEM) solver.  That will be a nontrivial task and I wouldn't suggest going in that direction unless you are an expert numerical analyst and have extensive experience in writing FEM code.

@KIRAN SAJJAN What you have listed under BCs is not enough.  You will also need to account for the condition expressed as equation (3) in the paper that you have cited.

That, however, won' t help.  There is no symbolic solution to that system of PDEs; I wish there were, but there isn't. If there were, then there wouldn't have been a point in publishing the paper that you have cited; they could have just published the symbolic solution instead.

The best you can hope for is a numerical solution.  But for that you need to specify the numerical values of the parameters, such as lambda, epsilon, Pr, etc.  Even after you do that, Maple's numerical PDE solver won't be able to help you find a solution because it is limited to solving evolution equations, that is, PDEs that depend on time, while yours isn't.  Furthermore, the condition (3) noted above is not a PDE so it needs a special treatment.  You need a specialized solver to solve the PDE numerically, and that would be a somewhat advanced research project.

@KIRAN SAJJAN The original system of PDEs is nonlinear, so we don't expect a symbolic solution which is Maple's main strength.  You wrote that you tried solving the system in Matlab.  That may be an option, but I don't know enough about Matlab's PDE solver to be of help.  You may want to ask in a Matlab forum.

As to

    In matlab also I have tied but for only 0 to 1 I'm getting
    or -1 to 0 only getting.

I can't tell what that means and it's not clear to me whether you fully understand the problem's setting.  The domain of the PDEs is the three-dimensional annular region bounded by the planes z=±h and a < r < b in cylindrical coordinates.

@KIRAN SAJJAN In the paper that you have cited, the author expands the solution into infinite series and thus reduces the PDEs to a set of ODEs (11) and boundary conditions (12).  The graphs shown are those of the solutions of that system of ODEs.

You wrote that you know how to solve ODEs, so you should be able to continue from there.

1. The equations indicate that the unknowns u, v, w, p, theta are functions of r and z, but the plots show these as functions of z only.  That's not consistent.

2. The PDEs model the motion of a fluid in a cylinder.  The radius of the cylinder is not given.  Also you need the parameter values of Pr, lambda, and epsilon.

3. Boundary conditions are given at the top and bottom (the flat parts) of the cylinder, but boundary conditions on the lateral surface of the cylinder (the curved part) is missing.

@Dkunb Most of this worksheet's functions execute in a fraction of a second on my (very old) computer.  There is not much to measure there.  Do you have any specific function in mind?

The slowest operations are in the calls to display() which take a few seconds.  To time each display() call, do

st := time():
display(... whatever ...);
time() - st;

That will print the number of seconds (and fractions of seconds) elapsed in executing the display().

As to your question regarding arranging the plots into 2 rows and three columns, do

display(< p1, p2, p3; p4, p5, p6>);

where the p1, p2, etc., are plotting commands, such as odeplot(...).  The semicolon after p3 separates the rows.

PS: Don't forget to adjust my worksheet according to the note titled Correction posted on May 10.

@Scot Gould Save mod for polynomial algebra.  It's safer using irem with numbers, as in

plot(irem(round(x),2), x=0..5,  thickness=4);

The coefficnet i (the imaginary unit) is missing from the set of equations for C in the worksheet that I have provided.  To correct, change the line

de_C := seq(diff(Cvec,tau) =~ lambda*exp(-y(tau)) * B . Cvec);

to

de_C := seq(I*diff(Cvec,tau) =~ lambda*exp(-y(tau)) * B . Cvec);

Then change the plotting range by setting t_final = 2.

@Dkunb You have specified y(0) = 20 and py(0) = -5.  That's a good start.  But we have N+2 first order differential equations, so need N+2 conditions to solve.  What are the rest?  And what is  the range of tau?

It will help if you supply a new worksheet where you specify the numerical values of all the parameters, such as m_x, lambda, N, and so on, a full set of initial conditions, and the vector E.  We don't need the differential equations themselves since we can read them in the image that you have posted.

@Dkunb The set of ODEs that you have shown at the end of the image that you have posted involves m_x, m_y, E_m, E_n, lambda.  Are these variables or constants?  What do we know about them?

The independent variable is tau.  I assume that it is a real (not complex) variable.  What is tau's range? Is it the entire real line?

What is the value of N?

Your system of ODEs involves N+2 unknowns.  To obtain a numerical solution, you need to supply N+2 initial or boundary conditions.  What are they?

I suppose one may extract the answers to these questions by sifting through your worksheet but that would be a Herculean task.  It will help if you could remove your attempt to a solution, and just present the parameter values as asked above.   We already know the system of equations from the image that you have posted.

@Paras31 For whatever it's worth, here is that diagram in tikz.  It takes fewer lines to code and the result looks more professional.

\documentclass{article}
\usepackage{tikz}
\usetikzlibrary{calc}
\begin{document}

\tikzset{Arrow/.style={->, very thick, red}}
\begin{tikzpicture}[>=stealth]
    \pgfmathsetmacro{\R}{3}
    \pgfmathsetmacro{\r}{0.13*\R}
    \pgfmathsetmacro{\H}{0.5*\R}
    \pgfmathsetmacro{\a}{70}
    \pgfmathsetmacro{\d}{\r*cos(\a)}
    \draw[->] (-1.2*\R,0) -- (1.2*\R,0) node[right] {$\Re$};
    \draw[->] (0, -1.2*\R) -- (0,1.2*\R) node[right] {$\Im$};
    \draw[red, semithick] (0,0) circle (\R);
    \draw[red, semithick]
        (0, \H) +(-\a:\r) arc(-\a:180+\a:\r)
        -- ($(0,-\H) + (180-\a:\r)$) arc (180-\a:360+\a:\r) -- cycle;
    \draw[Arrow] (0,\R)    -- +(-0.01,0);
    \draw[Arrow] (0,\H+\r) -- +(-0.01,0);
    \draw[Arrow] (0,-\R)    -- +( 0.01,0);
    \draw[Arrow] (0,-\H-\r) -- +( 0.01,0);
    \draw[Arrow] ( \d,0) -- +(0, 0.01);
    \draw[Arrow] (-\d,0) -- +(0,-0.01);
    \fill[red] (0,\H) circle (0.1*\r) (0,-\H) circle (0.1*\r);
    \draw[>=latex, <-] (0,\H) +(40:0.1*\r) [out=40,in=180] to +(20:2*\r)
        node[right] {$z=i$};
    \draw[>=latex, <-] (0,-\H) +(-40:0.1*\r) [out=-40,in=180] to +(-20:2*\r)
        node[right] {$z=-i$};
\end{tikzpicture}

\end{document}