For some older discussions of this topic you could see here and here.

Do you want a surface that interpolates the data (ie. the surface should pass through the z-values of the given data points), or a smooth surface that approximates the data (without necessarily passing through each data point)?

Sometimes the choice will depend on on whether you anticipate that the data values have some error due to (say, physical) measurement.

Once you've imported the data you might have a look at the example worksheet with help topic,
    examples,Interpolation_and_Smoothing

That example worksheet has a smoothing example for 3D data, using the ScatterPlot3D command (similar to Carl's fine Answer, which also shows you how to import the data).

Alternatively, in Maple 2018 the new Interpolation package can be used to produce an interpolating object (call it F) from you data triples. This F can then be called with other numeric x and y values, where that F(x,y) computes an interpolated numeric z value. And so you can then use that F inside 3D plotting commands like plot3d . The example worksheet contains an example of that too. There are a variety of methods and options for the interpolation.

@Markiyan Hirnyk The time to evaluate to floating-point is not of greatest importance here.

I also added an exact simplifcation, for the original expression evaluated at each of the roots. I agree with Carl that that is important consideration in such examples, although I don't see it in his Answer.

Another aspect that can be of importance in related questions is whether small float imaginary components are mere artefacts due to roundoff error, etc. And so the application of evalc is a useful kind of simplification technique to know. So while not directly relevant for this example it's still worth mentioning.

The application of evalc can also be also useful for the verification step. Eg, for Carl's solution for the roots, and using his f definition and code,

CLans:=simplify~(eval~(Z, {solve(evalc~({Re,Im}(f(Z))), {x,y}, explicit)})):

simplify(convert~(evalc(eval~(f(z),z=~(CLans))),expln));
                            {0}

@Carl Love That does produce a nice surface, and there's no seam (because the suface doubles up when c>0). On my machine the horizontal gridlines are not evenly spaced when c>0, however.

So please forgive me if I offer an alternative here, where the gridlines are evenly spaced (but there is a very slightly darker horizontal gridline at z=0 when c>=0 -- the choice being either that or a small number of chinks visible upon manual rotation).

Explore(plots:-display(  
          plot3d([abs(sqrt(z^2+c))],
                 theta=-Pi..Pi,z=`if`(c>=0,0,sqrt(abs(c)))..sqrt(2),
                 coords=cylindrical,thickness=0),
          plot3d([(sqrt(max(0,z^2+c)))],
                 theta=-Pi..Pi,z=-sqrt(2)..`if`(c<=0,-sqrt(abs(c)),-5e-4),
                 coords=cylindrical,thickness=0),
          view= [-2..2,-2..2,-sqrt(2)..sqrt(2)]),
        c=-1.0..1.0);

With style=surface all these gridline issues vanish, naturally.

@stefanv It works in the command line interface but not in a Worksheet in the Standard GUI.

@nm 

What I'm trying to say is that a program that forces repeated copying of (immutable structures) sets or lists by repeatedly assigning into entries (using whatever syntax) is already a poor program.

And any syntax that makes coding that easier, in any language, is a poor language feature.

A syntax which makes a poor language feature "more natural to use" is no virtue at all.

@Carl Love Reaching for the eliminate command seems the natural thing to do here. But it looks as if there's something to fix, because the following returns the empty set instead of restrictions on the remaining name k.

That is, I am seeing this,

restart;

eqs:={-6*c+(3/2)*c^2-2*b-3*b*c+(3/2)*b^2-3*a*c+(k^2)*b-b+(3/2)*(a^2)+(k^2)*a-a=0,
      -2*b-3*b*c+3*(b^2)-6*a*c+2*(k^2)*c-2*c-9*a*b+3*(k^2)*b-3*b+6*(a^2)+4*(k^2)*a-4*a=0,
      (3/2)*(b^2)-3*a*c+(k^2)*c-c-9*a*b+3*(k^2)*b-3*b+9*(a^2)+6*(k^2)*a-a=0,
      -3*a*b+(k^2)*b-b+6*(a^2)+4*(k^2)*a-4*a=0, 
      (3/2)*(a^2)+(k^2)*a-a=0}:

eliminate(eqs, {a,b,c}):

           [{a = (1/3)*k^2-1/3, b = 0, c = 0}, {}]

But I would have expected something with a restriction on k in one part of the result, say something like this,

    [{a = 0, b = 0, c = 0}, {}],

    [{a = (1/3)*k^2-1/3, b = 0, c = 4}, { k^2 - 1 }]

The difference between Maple 16.02 and Maple 17.02 is that frac(Pi^22) is no longer being returned as (the same) unevaluated frac call. I completely missed that before.

Eg,

restart;            
kernelopts(version);

      Maple 16.02, X86 64 LINUX, Nov 18 2012, Build ID 788210

frac(Pi^22);

                             22
                      frac(Pi  )

versus,

restart;            
kernelopts(version);

        Maple 17.02, X86 64 LINUX, Sep 5 2013, Build ID 872941

frac(Pi^22);

                        22
                       Pi   - 86556004191

So the difference is Library side (rather than any change in the underlying kernel-extension mechanism, AFAIK). So with the current control-flow (and with frac(Pi^22) returning just Pi^22-86556004439 ) it's not clear at what point the frac computation could be intercepted and Digits adequately raised internally (see Carl's exposition) so that evalf(frac(Pi^22)) would return an answer accurate to inbound Digits.

@ecterrab I wrote a comment here which I've now deleted. I missed the fact that, in Maple 17.02 onwards, frac(Pi^22) is no longer returning unevaluated as it did in Maple 16.02.

I like vv's idea of approaching the problem via a parameterization.

The real and imaginary components of the original can be plotted, for some insight. It's harder for implicitplot if those are two surfaces which intersect along a curve at height z=0. And it's harder still if either of them only touches z=0 without crossing it.

Using fsolve can work, provided there is only one relevant solution p for each b (or vice versa), and if fsolve can always find the desired other real coordinate value (eg, when called with its complex option).

Here a a few more brief considerations:

impl3.mw

@Earl The earliest version of the Explore command saw it get so-called special evaluation rules. In particular the first formal parameter of the Explore:-ModuleApply procedure (ie. the expression being explored) is specified uneval.

That brings the shallow benefit that one does not need to use single right-quotes (aka uneval quotes) to prevent the first argument from being prematurely evaluated. I call it shallow because it also prevents assigned names in the explored expression (here, the equations in sys, the assigned names within them, the procedure fric, etc) from resolving to their values unless special examination and evaluation is done internally by the procedure.

In Maple 2018 the internals of Explore go to some effort to try and resolve the assigned names. It mostly works¹.

In Maple 2016 and older there are some alternatives. One alternative is to create a wrapping procedure (which can, when called, get the needed evaluation via scoping), and this is why the ?examples,Explore example worksheet has several examples with simple wrapping procedures. Another alternative is to use some suitably placed uneval quotes, and wrap with a call to eval, and I attach that below.

382483three_reservoir_prob_ac_2016.mw

Sorry for the long-winded explanation.

¹⁾ After 25 years working with Maple I've become convinced that it is not possible to implement (as interpreted, Library code) exactly the same kind of evaluation of procedure arguments as the kernel can effect. So such extra efforts can never be made to work entirely.

An implementation of ithprime could utilize a stored collection of (suitably placed) primes and their ordinates.

When called with argument n the procedure would find the largest stored ordinal less than or equal to n. And then it could start the `isprime` from the corresponding stored prime (if n were not a stored value).

Storing 1 million integer[8] values takes only about 8Mb in memory. 

An appliable module could load the precomputed values from a binary file, when being defined.

Perhaps a fun project.

@Carl Love The setting is saved as a preference automatically.

Why do you want to make the whole legend be just a single string?

Why cannot you use either of these? (Are you using an old version of Maple?)

legend = [omega[0]=omega0]

legend = [`ω`[0]=omega0]

@David Sycamore By "argument" Tom is referring to n when ithprime(n) is called.