@brian bovril I spent some time this weekend on developing a procedure which serves to extrapolate a square's boundary data into its interior.  You will find the details in this MaplePrimes post.

The construction there assumes the values of the desired function are given on all four edges of the square.  It does not handle the cases, as in this current thread, where the derivatives are specified on some edges.  I believe that the method I have described in the link noted above can be extended to handle derivative data as well but I haven't attempted that.

@brian bovril I can explain that construction, however it is quite ad hoc and not particularly illuminating.  Here is an alternative solution (of course, the question does not have a unique answer) which is both neater and easier to explain.

Consider the Cartesian coordinate system r−z (r is horizontal, z is vertical) and the associated polar coordinates ρ−θ.  The function f(ρ,θ) = cos θ defined on the square [0,1]×[0,1] is 1 along the bottom edge and 0 along the left edge. That's the same as the function σ₁(r,z) = r/sqrt(r^2+z^2) in the Cartesian coordinates.  Through translation and reflection, we see that the function σ₂(r,z) ≡ σ₁(1-r,1-z), defined on the same square, is 1 along the top edge and 0 along the right edge.

The function σ₃(r,z) ≡ σ₂(r^2,z^2) defined on the larger square [−1,1]×[−1,1] is symmetric with respect to reflection about the coordinate axes, and therefore ∂σ/∂z=0 on the r axis, and  ∂σ/∂r=0 on the z axis.

We observe that the restriction of σ₃(r,z) to the fourth quadrant satisfies the desired boundary conditions.  The only trouble is that the original question places the domain in the first quadrant.  We remedy that through a translation and arrive at the final solution σ(r,z) ≡ σ₃(r,z−1).  Here are the details in Maple:

@ecterrab An excellent job!  As always, I am impressed by your quick response to suggestions made in this forum, and your ever-increasing contributions to Maple.

@rlopez It depends on what one assumes on "r".  Here are two scenarios:

1. We have w = f(b, phi, r), and g(b, phi, r)=0.  The latter may be solved, in principle, for r and thus give r = G(b, phi).  Then w = f(b, phi, G(b,phi) is a function of two independent variables b and phi.
2. We have w = f(b, phi) and g(b,phi)=0.  Then latter may be solved, in principle, for phi and thus give phi = G(b).  Then w = f(b, G(b)) is a function of a single variable b.  Alternatively, we may solve g(b,phi)=0 for b and get b=H(phi).  Then w = f(H(phi), phi) is a function of the single variable phi.

In my calculation I followed the first scenario, while in your calculation you follow the second scenario.  The difference is that in the first scenario r is a variable, while in the second scenario it is a constant. As vv commented, the problem is not clearly formulated.

@ecterrab I am happy to hear that further work is being done on this and I look forward to seeing and using the resulting product.

Thank you everybody for your solutions, comments, and insights.  Now that the problem is solved, I would like to tell you what it is all about.

Consider the standard quadratic equation a*x^2 + b*x + c = 0.  Suppose that the coefficients a, b, c, are selected randomly from a normal (Gaussian) distribution of mean zero and standard deviation sigma. We wish to determine the probability of complex roots.

Any selection a, b, c, of the coefficients corresponds to point (a,b,c) in a three-dimensional Cartesian coordinates system.  The implicitly defined surface b^2 - 4*a*c = 0 divides the space into two subsets, one corresponding to real roots, the other to complex roots.  Here is what the surface looks like:
plots:-implicitplot3d(b^2 - 4*a*c = 0,
    a=-4..4, b=-4..4, c=-4..4,
    style=surface, numpoints=10000,
    orientation=[-65,80,0]);

Aside: Actually that's the surface of a cone, but I am not going to exploit that property here.

It is evident that the region corresponding to complex roots is the union of two subregions
    { (a,b,c) : a > b^2/(4*c), c > 0  }  and  { (a,b,c) : a < b^2/(4*c), c < 0  }.
The two subregions are congruent.  Therefore it suffices to compute the probability of a point lying in one of them and then double the result.

@Mariusz Iwaniuk Thanks for your answer.  Actually, motivated by your success, I played around a bit by changing the orders of the integrations, and thus obtained an exact (non-numeric!) answer which of course agrees with yours:

@DJJerome1976 As to "it never dawned on me that shadebetween( ) accepts variable limits", the documentation may not be without blame.  The help page says:

shadebetween(f1, g1, x=a..b, y=c..d);
where "a, b, c, d are real constants".  (my emphasis)

Furthermore, due to the lack of sufficient documentation, it is not clear why the second of the following two commands fails to "shade in between" in Maple 2018.1:
plots:-shadebetween(0, x^2 + y^2, x=0..1, y=0..x);
plots:-shadebetween(0, x^2 + y^2, y=0..x, x=0..1);

Kitonum's plot3d() variant succeeds in both cases.

@mmcdara You have correctly described the classical Poiseuille flow.  In that context, the radial velocity v is zero and the axial velocity u is parabolic.

The problem originally stated in this thread is not quite the Poiseuille flow—the pipe leaks along its lateral surface as indicated by the boundary condition v=p+a at r=1.  Consequently, the radial velocity is nonzero, and as far as I can see, there is no solution in terms of elementary functions.  A numerical solution should work.

@mmcdara I think the pbar in the original statement is meant to indicate the average pressure.  Thus, the average pressure at z=-1 is zero and at z=+1 is p_a.

I don't know whether these are sufficient boundary conditions to make a well-posed problem.

@Ronan the last command in my worksheet contains the construction
    seq(frame(t), t=0..100, 0.4)

This produces animation frames over the time t=0 to t=100, in steps of 0.4.  You may change those numbers as you wish.

As to the blog that you referred to, I find the interface non-intuitive and therefore am having difficulty navigating through it.  I searched there for your statement regarding the ratio of the angular momentum and kinetic energy but was unable to locate it.  Can you post a direct URL?  The statement should involve more than just the ratio, since that ratio is not a dimensionless quantity.

@vv You are right; it's not difficult to compute moments of inertia of the T-handle exactly through direct integration.  Following your idea, in the attached worksheet I calculate the moments of inertia while accounting for the overlap between the two cylinders.  As expected, the results are quite close to the approximate calculation where the overlap is ignored.

T-handle-flip-extra.mw

@Carl Love Your function Z is a prototypical infinitely differentiable function which is not analytic.  As you have noted, all derivatives of Z at x=0 are zero.  As a consequence, a naive attempt at expanding Z(x) in Maclaurin series leads to an infinite sum of zeros which certainly does not equal Z(x).  This is because Z is not analytic, and that ties in to your comment about complex numbers.

A function closely related to your Z is
    Y := x -> 2.252283621*piecewise(abs(x) < 1, exp(-1/(1-x^2)), 0);
which is also infinitely differentiable (but not analytic), and has compact support.  The purpose of the coefficient 2.252...  is to normalize the function so that the area under the graph is 1.

Define Y_ε(x) = Y(x/ε). The function Y_ε is called a mollifier. For any function F in L²(R), let's write F_ε ≡ Y_ε o F for the convolution of Y_ε and F.  It can be shown that:

• F_ε is infinitely differentiable;
• F_ε converges to F in L² as ε goes to zero.

Note that although a function F in L² is generally not differentiable, the mollified version, F_ε, always is, and therefore the derivative F ^'_ε is well-defined.  Additionally, if F is such that F ^'_ε converges to something in L², let's say G, then G is called the weak derivative of F.

Sobolev in Russia and Friedrichs in the US capitalized on this and developed what is now the standard theory of weak solutions and regularity of PDEs.

The Wikipedia page https://en.wikipedia.org/wiki/Mollifier has some details and historical remarks on mollifiers.

Saving Maple's graphics to PDF produces a useless result.  I have reported this to the Tech Support multiple times over the last five years or so, and each time I have received acknowledgment of the existence of the problem.  I am still waiting for a fix.

To save a graphics produced by Maple's plot(), I right-click on the plot, select Export and then one of the several choices of graphics formats.  The options PNG, GIF, JPEG, EPS work fine—they produce graphics files whose bounding boxes correspond to the extents of the image.  Saving to PDF misbehaves—it produces the equivalent of an 8.5'' x 11'' paper and inserts the graphics somewhere near the upper left corner.

I am absolutely at a loss to see the utility of that.  What in the world is the use of an  8.5'' x 11'' export? Shouldn't exporting to PDF produce a natural bounding box as exporting to other image formats do? 

@Kitonum No need for pointplot().  With the output=operator option to dsolve(), we may use spacecurve() as you did in your first example:

restart;
sol:=dsolve({
    diff(x(t),t)=-sin(t),
    diff(y(t),t)=cos(t),
    diff(z(t),t)=0.3, x(0)=1,y(0)=0,z(0)=0},
    numeric, output=operator):
plots:-spacecurve(eval([x(t),y(t),z(t)], sol),
    t=0..4*Pi, colorscheme = [blue, green]);