Yes. Your last paragraph describes what's happening.

It's not specific to the simplify command.

Why are you so surprised by this?

What would be your comprehensive delineation between circumstances where NULL inside an expression sequence should get flattened out and circumstances where it should not?

It doesn't seem like a very good idea to have an application that relies on the initial state of currentdir().

But cannot your application's startup code set currentdir to something relative to, say, kernelopts(':-homedir') ?

What sort of i/o do you intend on trying within the MaplePlayer? (There's a possibility that the MaplePlayer would have to look into a temp directory in the FileTools sense.)

Right, an Execution Group or Document Block can only have a sole Task region, which can have worksheet content embedded into it. So each call to Tabulate (or Explore, or InsertContent, or ImageTools:-Embed, etc) executed within the same Execution Group simply replaces the content in that sole Task Region.

Here are two approaches to getting what you've described.

The first way involves keeping a running Matrix that accumulates the various results. It simply replaces the embedded tabulation each time that additional results are added to the bottom of that Matrix. This is easy if the tabulations all have the same number of columns.  Tab1.mw

The second way involves keeping a running copy of the XML representations of separate Tabulate calls, stored separately as, say, entries of a Vector. At each iteration it embeds all the tabulations accumulated so far (in order) as separate Tables. This is more suitable if the various tabulations are quite different from each other, eg. different number of columns. Tab2.mw 

The Tabulate command is really just a convenient way to build -- and, normally -- embed the XML representing a GUI Table. The undocumented option output=XML allows the Tabulate command to return the XML for the constructed Table without actually embedding it. So the second approach simply stores them, allowing you to replace the sole Task region for an Execution Group or Document Block with whatever combination of such previously constructed Tables as you see fit.

Hopefully you can adjust one of the two approaches to your purpose. I used Maple 2018.0.

You already had them combined by using plots[display], but their curves are so close that they overlap and you can only see the one on top.

One alternative is to plot one of them as a line and the other as a point-plot.

I increased h for Runge while plotting List as a point-plot, since otherwise the point symbols overlap and don't look so good. But you could also use the smaller h and just plot every 2nd (or 10th, or whatever) entry of List generated using the smaller h.

I also added a legend entry to each.

Another alternative is to plot them side by side.

question_ac.mws

It'd be helpful if you told us what exactly you wanted the final combined plot to look like.

I often find that it's more convenient to use so-called 2-argument eval than to use assign on the fsolve result, because the former allows you to immediately do further symbolic computations with the original variables.

But if you do the assign then you'd have to call unassign if you wanted to do further symbolic work with the original variable names.

The appeal of using assign is that formulas involving the names will get the numeric values directly.

Which approach you choose is up to you. Neither is incorrect. Which one is more convenient for you depends on what kind of further work you intend on doing.

restart;
eqs := [c + d = 520.9447091, c - b = 0.04148708484,
        d*a + c + b = 0.04147082666, b - c - d*a = -0.04148773634]:

res := fsolve(eqs);

                                                    -5                                      -8
     res := {c = 0.04147863000, b = -0.8454840000 10  , d = 520.9032305, a = 0.1250712151 10  }

assign(res);

sin(d) + c/b;

                     -4906.469192

cos(d) + b/a;

                     -6759.195990

unassign('a','b','c','d');

# Now you can do further work with the names as symbols.
fsolve( sin(d) = d - 0.5 );

                      1.497300389

# But now the formula doesn't immediately produce a numeric result.
cos(d) + b/a;

                     cos(d) + b/a

eval( cos(d) + b/a, res);

                     -6759.195990

And here is an alternate approach.

restart;

eqs := [c + d = 520.9447091, c - b = 0.04148708484,
        d*a + c + b = 0.04147082666, b - c - d*a = -0.04148773634]:

res := fsolve(eqs);

                                -8                                                          -5
     res := {a = 0.1250712151 10  , c = 0.04147863000, d = 520.9032305, b = -0.8454840000 10  }

eval( sin(d) + c/b, res );

                     -4906.469192

eval( cos(d) + b/a, res );

                     -6759.195990

fsolve( sin(d) = d - 0.5 );

                      1.497300389

cos(d) + b/a;

                      cos(d) + b/a

eval( cos(d) + b/a, res);

                      -6759.195990

You can approach this using operators and D, or expressions and diff. But note that neither will produce ten correct significant digits in the answer if you just use the default working precision of Digits=10.

Increasing the working precision can produce the desired result. But it's not any kind of proof.

restart;

f := 2*sin(x^2/2)-sin(1*x/2)^2:
df := diff(f,x):
dff := diff(f,x,x):

Digits := 20;
                          Digits := 20

r := fsolve(df, x=1..2);
                   r := 1.6870521236466276364

a := eval(dff, x=r);
                  a := -5.2779215898752826466

ans := evalf[10](a);
                      ans := -5.277921590

That may be adequate for you. Using Maple to try and prove correctness of an answer to this question is an interesting task, though it might be out of scope for you.

Two interesting aspects are: how accurate does r have to be such that substitution into the second derivative can robustly produce a result correct to ten significant digits, and what working precision might be needed for that substitution.

Even considering alone the second of those aspects is not basic.

restart;
f := 2*sin(x^2/2)-sin(1*x/2)^2:
df := diff(f,x):
dff := diff(f,x,x):
Digits := 40;
                          Digits := 40

r := fsolve(df, x=1..2);
         r := 1.687052123646627636437294496471862263832

for d from 10 to 20 do
  forget(evalf);
  Digits := d;
  s := shake(subs(x=r, dff));
  a := evalf[2*d]((op([1,2],s) + op([1,1],s))/2);
  t := evalf[2*d](op([1,2],s) - op([1,1],s));
  g := length(trunc(abs(a))) + floor(abs((log10(frac(abs(t))))));
  print(evalf[10](a), t);
  if g >= 10 then
    ans := evalf[10](a);
    break;
  end if;
end do:
                                          -8
                    -5.277921592, 2.255 10  
                                          -9
                    -5.277921589, 2.256 10  
                                         -10
                   -5.277921590, 2.257 10   

ans;
                          -5.277921590

Another way to consider the aspects I mentioned might be to successively compute r accurate to d digits, and then construct an INTERVAL(r-10^(-d), r+10^(-d)) and examine what could be done with that.

Note that simply using default Digits=10 working precision need not produce the answer correct to ten significant digits, whether done using operators or expressions. But both approaches could be used to find the answer correct to ten significant digits.

# Kitonum's suggestion
restart;
f := x->2*sin(x^2/2)-sin(x/2)^2:
x0:=fsolve(D(f)(x)=0, x=1..2);
                   x0 := 1.687052124
(D@@2)(f)(x0);
                      -5.277921591

restart;
f := 2*sin(x^2/2)-sin(x/2)^2:
r:=fsolve(diff(f,x), x=1..2);
                    r := 1.687052124
eval(diff(f,x,x), x=r);
                      -5.277921591

restart;
Digits:=20:
f := x->2*sin(x^2/2)-sin(x/2)^2:
x0:=fsolve(D(f)(x)=0, x=1..2);
                x0 := 1.6870521236466276364
(D@@2)(f)(x0);
                  -5.2779215898752826466
evalf[10](%);
                      -5.277921590
restart;
Digits:=20:
f := 2*sin(x^2/2)-sin(x/2)^2:
r:=fsolve(diff(f,x), x=1..2);
                r := 1.6870521236466276364
eval(diff(f,x,x), x=r);
                  -5.2779215898752826466
evalf[10](%);
                      -5.277921590

One way to do it is by using the plots:-densityplot command. But that can make the Maple GUI behave badly if the grid size is too large. And that allows you to shade the hue by the argument, but doesn't provide any nice mechanism for controlling the intensity or saturation layers.

Instead you can create a color image directly as a 3-dimensional Array.

The following can be made faster, but here I'm focusing on the process.

If you want to get fancy you could try and implement a piecewise structured adjustment to the intensity or saturation layers of the Array (as representing HSV format).

If you increase the grid size (m and n) a great deal then you'll be better off visualizing it with ImageTools:-Embed than with plot and its background option.

argument_plot.mw

See also this very old post. But ignore the followup Comments where I added plot axes, since that is done better and more easily in modern Maple by using the background option of the plot command (as I did above).

See also the MathApp accessible in the Help system under the topic ArgumentShading .

For comparison's sake, here is plots:-densityplot on the same example,

density_arg.mw

I could also mention that the color shading data can be pulled out of the plot structure returned by plots:-densityplot, as discussed in this old post.

The keystroke pair Ctl-t switches from math entry mode to text entry mode. That text will appear in black by default.

From text entry mode you can use Ctl-r to switch to active math entry mode, or you can use Shift-F5 to switch to inactive math entry mode.

You can switch back and forth, in order to produce 2D typeset math appearing inlined within a paragraph of text. The active math will get executed (when you hit the !!! icon on the menubar, say) while the inactive math won't get executed.

What you saw in red was 1D Maple Notation (plaintext) math, and not text.

You could put the values in a Matrix and then export using commands like ExportMatrix or ExcelTools:-Export . No need to copy & paste.

Dummy_5_ac.mw

If you put all the values in a list then you can use the select command to separate out the purely real values. Then you could pick off one of them (the first, or last, possibly after sorting them).

[edit] The first argument passed to select could also test for inclusion within a specific range. Basically, you make the first argument passed to select be a predicate which returns true or false according to some condition of your choice. The third and additional arguments get passed as options to the first argument. See also Kitonum's answer below.

note: The fsolve example below is deliberately artifical. It's just an example. (If you were really using fsolve then you could specify a real range, or drop the complex option, etc.

For example,

restart;

f := a -> fsolve(x^3+a*x+1,x,complex):

f(2.3);

   -0.405741100177881, 0.202870550088940 - 1.55674962029228 I,
   0.202870550088940 + 1.55674962029228 I

sel := a -> select(type,[f(a)],realcons)[1]:

sel(2.3);

               -0.405741100177881

V := Vector(7):
for i from 1 to 7 do
  V[i] := sel(10*sin(i));
end do:

V;

                [-0.118641053844103]
                [                  ]
                [-0.109829320596285]
                [                  ]
                [-0.574354935221908]
                [                  ]
                [-2.81483317335910 ]
                [                  ]
                [-3.14753110820927 ]
                [                  ]
                [-1.82790314819886 ]
                [                  ]
                [-0.151678953479572]

Vector(7, i->sel(10*sin(i)));

                [-0.118641053844103]
                [                  ]
                [-0.109829320596285]
                [                  ]
                [-0.574354935221908]
                [                  ]
                [-2.81483317335910 ]
                [                  ]
                [-3.14753110820927 ]
                [                  ]
                [-1.82790314819886 ]
                [                  ]
                [-0.151678953479572]

For a univariate polynomial (after substituion for B and k) the sorting done under evalf(RootOf(...)) may well interfere with the notion that any single indexed RootOf represents a particular root (whatever than means). So there may not be much additional conceptual harm in allowing fsolve to compute them all and sort them as it does.

This is pretty quick.

Download file2_acer.mw

And, with option discont=true 

Download file2_acer_discont.mw

You could regain quite a bit of the original's speed by putting the discont=true option only plt3.  file2_acer_discont_plt3.mw

Suppose that your Matrix has already been assigned to the name M.

Create a new Matrix where entries of M which equal zero are replaced by undefined.

Then call plots:-matrixplot or plots:-surfata on that new Matrix.

newM := map(u->`if`(u=0,undefined,u), M):

plots:-matrixplot(newM);

plots:-surfdata(newM);

I'm not really sure what your questions are. Are you looking to get something like this?

simplify_question_ac.mws

It's straightforward to construct re-usable procedures that map CDF and Quantile over z as a list (or Vector, etc), as well as handle scalar z.

I have reversed Maple's usual meaning for the numeric::truefalse option here. So rational inputs are treated with the floating-point method (and exact symbolic handling is allowed by specifying numeric=false).

This isn't restricted to the Normal distribution. You could make X a RandomVariable of some other distribution.

Download cdf_map.mw

[edit] I could also mention that you can also use the elementwise ~ on the CDF and Quantile commands. (One reason I constructed the above procedures is so that the numeric option behavior was reversed, to use floats by default because that might be more Matlab-ish.) You might find this just as easy,

restart;
with(Statistics):

X := RandomVariable(Normal(0,1)):

CDF~(X, [2.7102, 2.2, 1.7]);
   [0.996637867609054, 0.986096552486501, 0.955434537241457]

Quantile~(X, %);
     [2.71020000000130, 2.19999999999921, 1.70000000000084]

CDF(X, 2.7102);
                       0.996637867609054

CDF~(X, 2.7102);
                       0.996637867609054

restart;
with(Statistics):

CDF~(RandomVariable(Normal(0,1)), [2.7102, 2.2, 1.7]);

   [0.996637867609054, 0.986096552486501, 0.955434537241457]

Quantile~(RandomVariable(Normal(0,1)), %);

     [2.71020000000130, 2.19999999999921, 1.70000000000084]

It seems that the solution obtained at lower values of n can provide a useful initial point for a higher value of n.

Below it's done in two ways.

The first uses this bootstrapping approach for values of n=7,13,17,18,19,20.

The second does it in a loop, using all n from 2 to 20. (slower overall, but perhaps more robust!?)

Download solveSYS_ac.mw

If we plot the values for each variable, in the solution as n changes, we might observe that all the curves appear linear except for p[j,fail] , j=1..4, which are monotomic and could be interpolated. I would guess that reasonably close initial points could be generated for the higher n, based on a smaller number of the solutions for lower n (which compute more quickly).
solveSYS_ac2.mw
I would not be surprised if symbolic examination of the equations could lead to a non-numeric or hybrid symbolic-numeric approach where the solutions (dependent on n) were constructed very quickly. I don't have time for that.