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.

Download nz.mw

It turns out that increasing Digits up front may not be necessary here. But I include it anyway, since it may give you an idea in the future.

The help page for parse describes that the offset option can be supplied alongside the statement option, and that the lastread option allows for a running update to a specified name to hold the position least read.

Doing that in a loop allows for multiple statements to be parsed from a single string. Or you can construct a procedure to execute such a loop.

That is more robust than simply splitting the string at colons and/or semicolons, since those could appear in locations where the parsing should not be split.

An extra choice is to try and discern (and return) information about what (previously unassigned) names might have been newly assigned by this process. I realize that this may be of less relevance to your purpose, since your script that constructs the strings likely already knows to what names it will assign.

Download parsestatements.mw

By "newly assigned" I mean unassigned before the call. I don't mean re-assigned.

The error message describes the situation: the Int call is not supported by evalhf.

Did you just want to make it faster? Note that the quadrature scheme can utilize evalhf for the evaluations of the integrand.

Download evalhf_int_ac.mw

Or do you have other computations to go inside procedure f1, which you also want done under evalhf?

The subs command does syntactic substituion without evaluation (unless the command name is indexed by the key eval).

The eval command does a more mathematical substitution, in the sense that it has additional knowledge of the mathematics germane to some substitutions.

Often the two produce the same result. Sometimes one specifically wants to avoid evaluation prior to being able to do subsequent manipulations of the result.

Sometimes subs is a bit faster. Relying on that only makes sense up to the situations where one doesn't need the mathematical knowledge.

subs(x=1.2, sin(x));
                         sin(1.2)

subs[eval](x=1.2, sin(x));
                       0.9320390860

eval(sin(x), x=1.2);
                       0.9320390860

There are also differences in how substitutions are done. The subs command allows you to do them together or in sequence. Eg,

subs([a=b,b=c], [a,b]);
                          [b, c]

subs(a=b,b=c, [a,b]);
                          [c, c]

subs(b=c,a=b, [a,b]);
                          [b, c]

eval([a,b], [a=b,b=c]);
                          [b, c]

I'm not sure which aspects are most relevant to your query. There are other subtleties or examples which might interest you.