@mmcdara Thanks for the kind words.

The calculation in my original worksheet extends almost without change to the general case of rotating a any curve (not necessarily a straight line) about a given line.  This is implemented as a procedure in the attached worksheet. 

Here is a sample obtained by rotating the helix <cos(t), sin(t), t/4> (shown in green) about the line <t, 0, t>. (shown in red). I have removed a part of the surface to make it possible to see the inside.

Download worksheet revolve-curve.mw

Although you have a system of three equations, you seem to be interested in plotting the system's vector field in the [r,phi] plane.  That is not possible in general. The vector field is 3-dimensional, so even when restricted to the the [r,phi] plane, the vectors will point outside that plane in general.  It's not a 2-dimensional picture.

@Vortex Unfortunately I know nothing about this area of physics.  I can only surmise from what you have said that f is a function of frequency, and the integration limits a and b are lower and upper bounds on the frequencies.  I have no idea what q represents.

It seems to me that the complete resolution of the question that you are asking depends on a full understanding of the physical phenomena that it models.  I am afraid I don't have that understanding.

@Vortex I am gathering from your response that the argument q ranges over the interval [a,b].

As p and q range over the interval [a,b],  the argument of f(p) ranges over [a,b], while the argument of f(p+q) ranges over [2*a,2*b].  So the domain of f is the union of those intervals.  For instance, if a=10, b=12, then your equation demands that f be defined over the interval [10,12] as well as the interval [20,24].

Under that scenario, it is reasonable to expect that one prescribes the values of f on [20,24] and seeks the values of f over [10,12].  Or perhaps one prescribes the values of f over [10,12] and seeks the values of f over [20,24].  The choice depends on the application in which this equation arises.

Here I have chosen a and b on purpose so that the intervals [a,b] and [2*a,2*b] do not overlap.  If they do, as in the case a=1, b=100 that you have cited as an example, the situation is more complex and the assumptions on where f is prescribed and where it needs to be calculated is not so clear.  Again, the source of the application should clarify the choice.

1. We see that f is not uniquely defined because if f is a solution, then so is f+C, where C is any constant.
2. I have no idea how to go about finding a symbolic solution.  My first attempt for finding a numerical solution would be to discretize the intervals (a,b) and (a+q,b+q), and look for a piecewise constant (or perhaps piecewise linear) f on those intervals.  The details will depend on how the size of q relates to a and b.  Any information about numerical values of the parameters?  Any information about the sign of f?

@rockyicer To remove the RootOf expressions, you may apply allvalues() to V1, as in

V1_alt := allvalues(V1);

That produces two solutions which you can examine through V1_alt[1] and V1_alt[2].  Neither matches your hand calculations.  I didn't have the patience to tracked down the source of discrepancy.

Consider assigning numerical values to all the parameters, evalaute your solution versus Maple's, and see whether they are actually the same thing in disguise.  If they aren't, then evaluate the set of your five equations by substituting your solution and Maple's.  See which one actually satisfies the equations.

@Carl Love Thanks for your quick reply.   I hope that someone from Maplesoft reads this and adds the caveat to the help page.

Furthermore, the help page should be explicit about the meanings of the spherical phi and theta angles.  Are they the physics phi and theta or the math phi and theta? I can figure out the answer with some experimentation but that's not ideal.

It will help if you upload the worksheet that produced those graphs. 

@Craber720  In my instructions I wrote:
    ...adapt it to the case of a cubic (3rd degree) polynomial that goes through four prescribed points. 

You have prescribed many more than just four points.  You cannot play fast and loose like that in mathematics.

Carl Love offered several interpretations of your vaguely stated question.  You need to be much more specific about what you need in order to get a concrete answer.  If this is a homework problem, don't muddle it by paraphrasing.  State what the original question asks.

@Oliveira As mmcdara has pointed out, this problem is beyond the capabilities of Maple's numerical PDE solver.  Maple's symbolic solver, however, should be able to produce a series solution.  Unfortunately it doesn't due to what appears to be an bug somewhere in the implementation.  This streamlined version of your PDE shows the issue.

PS: Note the minus signs on the right-hand sides of bc2 and bc3.  Those minuses are missing in your worksheet.

Download heat-conduction-in-cylinder.mw

That equation on its own makes not much sense.  At best it can be interpreted as a ordinary differential equation in the unknown x(y), where f(x,y) is known.

To get a useful feedback, you can begin by telling us where that equation comes from.

@tomleslie you wrote:
    The first thing to realise is that pdsolve(...., numeric) does not return any information
    about the derivatives of the dependent function T(y,t).

To get pdsolve() to return information about the y derivative of T(y,t), introduce that derivative as an unknown. Thus, replace all occurrences of diff(T(y,t),y) in PDE1 with S(y,t), and arrive at the system of two PDEs in the unknowns T(y,t) and S(y,t), like this:

PDE1 := Pr*(diff(T(y, t), t)-Ree*S(y,t)) = (1+Nr*(T(y, t)+1)^3)*diff(S(y, t), y)+3*Nr*(T(y, t)+1)^2*S(y,t)^2;
PDE2 := diff(T(y, t), y) = S(y,t);
ICandBC := {T(1, t) = 1, T(y, 0) = 1, S(0, t) = T(0, t)};

Then in the OP's original code (or your code) replace PDE[i] := {PDE1} with PDE[i] := {PDE1,PDE2} and now you can plot S(y,t) as you wish.

@gkokovidis That's strange because the symbol is available in both Document and Worksheet modes in my Maple 2022.2.  Can it be it is OS-dependent?  Mine is Linux.

1. You wrote:
   The conversion in maple to convert an ".mw" file in xmaple for a text file to be executed with maple command does not work.
   I have no problem exporting an *.mw file from within Maple as a plain text *.mpl file.  Perhaps I have misunderstood you.
    
2. Maple's *.mw files come in at least two different flavors, corresponding to Worksheet or Document formats, or more precisely, whether the input was performed as 1D or 2D text.  If the input is 1D, then it should be possible, with some effort, to convert the *.mw file to a plain text *.mpl file with the help of external tools such as sed/awk/perl. If the input is 2D, you will need Maple for conversion.

Although Maple's PDE solver is not equipped to solve your system of PDEs, that system is well-posed and can be solved both through the method of finite differences and the method of lines.  Depending on how familiar you are with such methods and Maple programming, the implementation may be a quick job or a long project.  As a starting point, introduce a new variable, u, defined through

u(xi,tau) = v(xi,tau) - phi__7*diff(v(xi,tau), xi,xi);

then use that to eliminate v from the equations.