@imparter Upload your worksheet so that I may modify it.

@imparter I don't see a question in what you have written, so I can't tell what sort of response will be of help to you.  But for one thing, I see that you have not included the linestyle=dash in your code that I had suggested.  Did you try that?

@imparter The plot option linestyle=dash will plot with dashed lines.  Additionally, you may specify a list of colors for the graphs if you don't like the default colors.  So you may replace your current plot commands with

my_plot_1 := plot([%], r=0..h(z), color=["Red", "Green", "Blue", "Cyan"], linestyle=solid);
...
my_plot_2 := plot([%], r=0..h(z), color=["Red", "Green", "Blue", "Cyan"], linestyle=dash);

and then superimpose the two plots with

plots:-display(my_plot_1, my_plot_2);

Posting a picture of your worksheet is useless because I can look at it but cannot execute it. To get helpful responses, upload your Maple worksheet.  How?  In the message editor's toolbar, click on the big fat green arrow to upload your worksheet.

@imparter I don't know much about HPM and honestly, am not much interested in it.  From the page of calculations that you have supplied, it seems to me that the calculations are elementary and can be done easily by hand, if you know what you are doing, which I don't.  You don't need Maple for those calculations.

@imparter The purpose of your calculations is not clear to me.

1. You are looking for an approximate solution.  It seems to me that you are seeking an asymptotic expansion in terms of the parameter p.  That will be good for small values of p. Is that what you are after?
2. You say (r^2 - h^2)/4 is an approximate solution.  That doesn't look correct to me.  If you plug that expression into the differential equaiton, the result is 1/2.  It should be something like zero.  Do you have a typo there?
3. You have not specified boundary conditions.  Are you looking for a general solution?  You have defined an expression labeled ibvc which looks like an attempt to specify boundary conditions, but that doesn't do that.

@emendes In the file charsets/v2.0/charsets we have 100+ occurrences of the double-quote character (").  In modern Maple this has been replaced with the percent sign (%).  So the first thing to do is to replace all " with % in the file v2.0/charsets.  (Use your editor's global find/replace for that.)

There is one additional error to fix.  At the very end of v2.0/charsets we have

save `charsets.m`;

That's not good.  The syntax for Maple's save command is

save name1, name2, ..., namek, fileName

The variable names name1, etc., are missing.  I don't know what this application is supposed to do, so I cannot supply those variable names.  Maybe you can.  This may be as simple as uncommenting the commented-out save command near the end of v2.0/charsets.  Try it.  I haven't.

You have w(r,z) for the unknown of the differential equation.  However z appears neither in the differential equation nor in the boundary conditions.  Shouldn't that be just w(r)?

1. I know nothing about the Ritt algorithm, however Maple 2024 (perhaps earlier too) has a DifferentialPrem function which seems to implement that algorithm.  Type "?Ritt" to go to the help page.

2. If that doesn't meet your needs and you still want to do it through Wang's package, it would be a good idea to let people know where you found that package so that they can see if it is still usable.

@imparter Here is the corrected worksheet.  I have replaced all occurrences of h by h(z).  I have also changed the Maple code that defines h(z) into a Maple proc because it works cleaner that way.

Aside: I see that you tend to terminate most of your statements with colons instead of a semicolons. Colons suppress Maple's output.  I would be very uncomfortable doing that.  I want to see the output of individual commands in order to confirm that at every step Maple does what I expect it to do.

Error-fixed.mw

@imparter You have shown the expected graphs of w(r) at z=1.  To get those graphs, set z=1 in your formula for h, calculate h, then use that h to plot the graphs as I shown in my earlier message.

By the way, in my code I have "M in [2, 5, 7 ,9]".   That's not correct.  The numbers in the square brackets are expected to be the values of M^2, not M.  To correct that, replace "M in [2, 5, 7 ,9]" with "M in sqrt~([2, 5, 7 ,9])".

@imparter I do realize that h is not a constant.  That's why I presented two samples with h=1 and h=0.8.  You may pick any h that you wish.  Or if you prefer, pick a value of z first, and then use your formula to calculate the corresponding value of h, and then use that h to plot the graphs.

I don't understand what you mean by "how to combine all those in one graphs for different range".  What function do you want to plot?

@imparter 

Domain is 0 to 1?  You must be kidding  :-(

Your unknowns w(r,z) and theta(r,z) are functions of two variables, r and z.  So your domain is a region in the two dimensional r,z plane.

To set up the PDEs properly, you need to

1. specify the geometric shape of that two-dimensional region;
2.  identify that region's boundary (precisely, not just hand-waving);
3.  specify the conditions that the unknwons w and theta take on that boundary.

@imparter You are attempting to solve a pair of PDEs in the unknowns w(r,z) and theta(r,z)..  You specify the boundary conditions as

conds:= w(h(z),0)=0, theta(h(z),0)=0, D[1](w)(r,0)=0, D[1](theta)(r,0)=0;

What is the domain of the functions w(r,z) and theta(r,z)?  You should specify the domain clearly and unambiguously before attempting to prescribe a meaningful set of boundary conditions.

1. The definition of h involves delta.  You need to assign a numerical value to delta.

2.  In the definition of pdes you have (DP+theta*x4*sin(alpha)).  Shouldn't that theta be theta(r,z)?