@Ioannis This worksheet adds labels.

Download: paraboloids-with-labels.mw

@tomleslie The original question asks for "barrel" shapes with circular profiles of radius R, consistent with the equation r = a - R*(1 - sin(phi)).  Your construction produces barrels with sinusoid shapes which is not the same thing.  The OP should refer to Kitonum's construction which produces proper circular shapes.

@Ioannis Here are a few issues for you to consider.

1. The picture that you have supplied in your original post does not agree with the functions that you have included there.  Here is the picture that corresponds to your functions:You  need do decide if you want to change to functions to agree with the picture, or to change the picture to agree with the functions.
2. As you have observed, your arrows do not come out correctly. That's because you are plotting one of the arrows as a binormal vector.  But that's not what you want.  The binormal vector pertains to curves, not to surfaces.  The tangent vectors that you are looking for need to be calculated as the derivatives of S1 with respect to eta and phi.
   
   It is perhaps possible to calculate and plot those vectors through the commands available in the Student:-VectorCalculus or Student[MultivariateCalculus] packages but I cannot help you there since I don't use them.  I find it easier to do such things without the Student packages as I have done to plot the mage shown above.  I am not giving its details right now because I cannot tell whether that's what you are looking for.
3. As to: "Also in the picture included the line (φ=const)".  But there is not just one (φ=const) line; there are infinitely many of them.  All the radial lines in the image shown above are (φ=const) lines.

@Ioannis You need to be more specific about "the command display doesn't work".  Does it result in an error, or you don't like what it produces?

As to "In the picture we have two paraboloidal surfaces!!! This is only the left!", I disagree. In the picture that I included in my earlier reply, you can clearly see both paraboloids.  Aren't you getting the same things?

@Ioannis What don't you like about this?
 

plots:-display(LH,LG);

See https://www.mapleprimes.com/questions/231331-How-To-Plot-A-3D-Graph

@lcz what do you expect this sum to be:

sum(milk, i=1..5);

That sum is 5*milk.  That's because the i in milk is not an index.  The expression 'milk' represents an atomic symbol; it is not a composite.  In the same way x__i is an atomic symbol; it is not a word with the index i tacked on.  See the suggestions in kitonum's post on how to create indexed symbols.

I have not examined the intricacies of your specific example but as a general answer to your main question, you may consider this.

We have an ODE x' = f(x,t,p), where x=x(t) is an n-vector (n=2 in your case) and p is  a parameter.  Let x1 and x2 be the solutions to the ODE corresponding to the parameter values p1 and p2.  We want to plot the difference x1(t) - x2(t).

Rather than solving the ODE twice, I would solve the system of ODEs

x1' = f(x1, t, p1)
x2' = f(x2, t, p2)

simultaneously as a system of 2n equations in 2n unknowns.  Then plotting the difference should pose no problem.

What you have described is quite specialized and I certainly don't have the necessary knowledge to construct a mathematical model of it.  If you have constructed a mathematical model and your difficulty is in its implementation in Maple, then post your worksheet and point out where Maple help is needed.

You will need to be more specific.  Plotting this function over 0..1000 is instantaneous:

plot(sin(t/100), t=0..1000);

@nm In Linux, Maple is launched from a script file in Maple's installation directory.  The script file processes command-line options and sets some environment variables before launching the binary executable.

I must have learned about the "–standalone" option from browsing that script file because the option does not seem to be documented anywhere.  Perhaps that's because it's specific to Linux?  I don't know.  There are some other useful undocumented options there.  The "–j" option, for instance, may be used to increase the Java heap memory allocation size which is sometimes necessary when doing extensive graphics.

I don't know about Windows, but in Linux the command-line option "–standalone" starts an independent Maple process.  There ought to be something equivalent in Windows.

I know how to solve that problem but I would rather not do your homework for you.

I will be happy to show you how to do it in Maple if you show me how you would go about solving it by hand.

@Mo_Jalal To plot u(x,1) we do:
pdsol:-plot(u, t=1);

To obtain solution values, we define U:
U := pdsol:-value(u);
Then U(0.8, 0.2) returns
            [x = 0.8, t = 0.2, u = 0.161382120212466312]
There are a lot of other ways.  Have a look at the help page on pdsolve,numerc.

As you have suspected, those initial conditions are totally wrong.  To pose a proper set of initial and boundary conditions, you need to provide more information about the problem that these equations are supposed to model.  Here I will provide one set of plausible assumptions but these may be entirely different from what you really need.

The main trick here is the introduction of the new variables u_t, u_tt, etc., which converts your original 3rd order system to a first order system.

Download mw.mw