restart;

X := Statistics:-RandomVariable(Normal(0, 1)):
pp:=[attributes(X)][3]:
pp:-ParentName;
                       NormalDistribution
pp:-Parameters;
                             [0, 1]

Y := Statistics:-RandomVariable(Normal(5, 17)):
pp:=[attributes(Y)][3]:
pp:-ParentName;
                       NormalDistribution
pp:-Parameters;
                            [5, 17]

Z := Statistics:-RandomVariable(Exponential(3)):
pp:=[attributes(Z)][3]:
pp:-ParentName;
                    ExponentialDistribution
pp:-Parameters;
                              [3]

You can easily construct your own float[8] Array to denote the color, using either HUE or RGB specification.

Download color_pointplot.mw

On MS-Windows I add the following undocumented option (following a blank space) to the launch command in the Properties of the desktop icon,

  -standalone

The main effect, apparently, is that the icon subsequently launches wholly separate java virtual machines.

If one of these GUI instances crashes (while printing a huge result, rotating a huge plot, waiting on a crashed Maple kernel, etc) then the other GUI instance remains responsive and its worksheets can still be saved properly and promptly, etc.

Since I quite often run multiple GUI sessions, often with examples that can crash the GUI, this has been very useful in saving my other work in progress.

My machine has adequate RAM, and I don't care that separate GUI JVMs can take up 300MB footprint each while running. It is much more valuable to me that I don't lose work in progress.

(Sorry, I did not notice this Question before now.)

One way to disable the discontinuity check done by evalf(Int(...)) is to use operator form for the integrand. And fsolve prefers everything in operator form when it detects such (even if not really necessary for a numeric scalar like 0.5 on the rhs).

restart;
fe := Int(x -> sqrt(1+(-5.557990765*sin(5.557990765*x)-7.3*cos(5.557990765*x)
               -5.6*sinh(5.557990765*x)+7.3*cosh(5.557990765*x))^2), 0 .. al) = 0.5:

fsolve(unapply((lhs-rhs)(fe), al), 0..1);

                          0.1550088024

That works even when the original int call had already been formed for the given example: We can recast the int call with integrand in expression form to an Int call with integrand in operator form.

restart;
fe := int(sqrt(1+(-5.557990765*sin(5.557990765*x)-7.3*cos(5.557990765*x)
          -5.6*sinh(5.557990765*x)+7.3*cosh(5.557990765*x))^2), x=0 .. al) = 0.5:

F := subsindets(fe, specfunc({int,Int}),
                u->Int(unapply(op(1,u),lhs(op(2,u))),rhs(op(2,u)),op(3..,u))):

fsolve(unapply((lhs-rhs)(F), al), 0..1);

                          0.1550088024

Now some comments about Maple's numeric integration.

It is long-standing and unfortunate behaviour that evalf(Int(...)) (ie, numeric integration) offers no all-round easy option to disable its examination of the integrand for discontinuities.

The important 1D numeric quadrature schemes utilized have error estimates that rely on certain degrees of smoothness, and thus they should often split the domain in the presence of a finite number of discontinuities. But there are lots of examples where the check for discontinuity consumes enormous resources, and where that check is not necessary or wanted (because there are no such points).

The typical unfortunate scenario occurs when a symbolic discontinutiy check on an expression containing float coefficients causes a symbolic expansion which runs amok when the floats get converted to exact rationals with large numerator or denominator.

Sometimes the check can be avoided by specifiying the method=_d01ajc for evalf(Int(...)) example. But that method requires Digits<=15, and when the calling environment has default Digits=10 then fsolve may try and resolve the question of convergence (or a root candidate) at yet higher working precision.

While it would indeed be useful if fsolve were to accept an accuracy tolerance option, the runaway crash in this particular example is due (IMO) to a bug in evalf(Int(...)) and not in fsolve.

In my opinion it is a very longstanding and serious bug that evalf(Int(...)) offers no option to forcibly disable its discontinuity test. I have answered dozens of Mapleprimes questions over the years where evalf(Int(...)) showed serious performance problems for users, where it was not originally obvious that the discontinuity check was the culprit.

I have seen several examples that don't involve fsolve at all, where the evalf(Int(...)) misbehaviour causes unnecessary very slow performance or even a crash. Some of these are plotting examples. Some are simply single instances of a numeric integration. The central problem is not in fsolve.

Consider this example, which runs amok (don't try it if you don't want a crash):

restart;
Q := Int(sqrt(1+(-5.557990765*sin(5.557990765*x)-7.3*cos(5.557990765*x)
         -5.6*sinh(5.557990765*x)+7.3*cosh(5.557990765*x))^2), x=0 .. al):
evalf[25](eval(Q,al=0.3));

The way that I disable the discontinuity check at the start of this Answer is to use operator form for the integrand. That can become very tricky when the problem involves nested integrals or integrands that depend on (scoped) parameters, etc. It can become so tricky that it really cannot be considered a practical alternative to a (currently missing) dedicated option to prevent the discontinuity check.

Here is the first plot obtained in the worksheet below, but without the blue plane. It has a smooth lower edge and no missing portion (for larger y values).

An upper bound y=20 was supplied. For the first approach below you could increase that (and the `view`, accordingly). An upper bound for the z was also supplied, in the `view` option, which you could also adjust to taste.

Download plot3d_inequality.mw

Here are a few ways, each involving just a few calls to LinearAlgebra:-GenerateMatrix.

I didn't keep exactly the same names as you have for Vector X, etc.

I include the situation where the final LHS is left in the form Matrix.diff(X) , as opposed to being fully explicit.

(But I would not be surprised if there were utilities to help in DEtools or PDEtools.)

[edit] I have included a way to print the Matrix equations (using either X.DX, or DX explicitly on the LHS) with inert multiplication. In case you were after such a visual display. For the case of the explicit LHS I've separated the inverse of X using its determinant factor and its adjoint.

Download various2.mw

I don't see why it's necessary to call expand here.

This produces a result close to your supplied target expression.

By the way, it looks like you gave a target expression with a wrong sign on one L1/n12 subterm.

Download collect_diff_ac.mw

You haven't told us what numeric value your "to" has.

So here are a couple of ways to do it purely as numeric quadrature, quickly, with the domain being specified by the parameter of a procedure.

Download evalf_int_example.mw

Don't use linalg, which is deprecated. You can extract any row of a Matrix with simple indexing.

You don't need to call print in order to display the scalar value assigned to a name.

I have used the parameter nr of aFshear to denote which row to use.

Download zip_assign.mw

You have not yet given numeric values to Nturns and Nsets.

The unapply command can make an operator from a previously computed expression (or a name to which such has been assigned).

The way you were trying to define M muddled the concepts of the formal parameters z and r of the operator with the names z and r in the expression previously computed and assigned to Bzsum. This is a common misunderstanding. In contrast, your simpler example (z,r)->z*r provides the body of the operator literally.

Download unapply_again.mw

As a programming excercise perhaps it's easier to start off with examples that you know have only rational roots, and then afterwards move on to more general randomly generated examples (whose roots may not all be rational).

In other words, start off by figuring out how to generate and test the set of candidates using examples for which you know a non-empty result should be found.

So, below I first deal with the case that all roots are rational. Then I move on to the case that possibly only one root is rational, to check the same technique. That seems like a reasonable beginning for you, before you try deal with the general case of randonly generated coefficients.

If you attack the general case up front then you may get muddled about whether no rational roots exist versus whether your code cannot find any properly.

Download ratroots.mw

I am guessing that your set may contain terms like a+b and that you would want that to be treated like the case of a constant or solitary parameter.

restart;

func := S -> not andmap(has,S,{x,y}):

func( {a*x^3-1, b*y+x, -1, c} );

                true

func( {a*x^3-1, b*y+x, c} );

                true

func( {a*x^3-1, b*y+x, -1} );

                 true

func( {a*x^3-1, b*y+x, b+a} );

                 true

func( {a*x^3-1, b*y+x} );

                 false

If you want to test only for constants or stand-alone parameters then it could be done by making the test check whether the entry s is constant or member(s,[a,b,c]) . Thus this treats the polynomial a+b differently from above. (I'm guessing that this is not what you want.)

params := [a,b,c]:

func3 := S -> ormap(s->s::constant
                       or member(s,params),
                    [op(S)]):

func( {a*x^3-1, b*y+x, -1, c} );

                true

func( {a*x^3-1, b*y+x, c} );

                true

func( {a*x^3-1, b*y+x, -1} );

                 true

func( {a*x^3-1, b*y+x} );

                 false

@EB1000 As of Maple 2018 the Library-generated commands appear only in the right side-panel (a.k.a. context-panel) and not in the popup.

The last item in the popup you've shown is "Open Context Panel for more...". You can click on that, or open the right side panel using the double side-chevron icon at right of the regular GUI menubars.

See the help page with Topic updates,Maple2018,ContextPanel

Did you obtain this Maple package from here? If so then you really ought to give proper attribution.

Do you have to have PHCpack installed somewhere on your machine?

Before you can run an example from within Maple do not you have to call the comand phc:-setPHCloc from the Maple package (which you reference in your Question here)? Perhaps you are supposed to call that with the location of PHCpack on you machine, perhaps at the start of each Maple session in which you intend on using this phc package.

I downloaded a precompiled 64bit Linux binary executable of PHCpack from here. I unpacked that to obtain an executable binary named phc . I put that in the same subdirectory as I use below to run the Maple script that creates the Library archive (and to which I save this Maple package).

The I ran the following in Maple 2017.2,

restart;

currentdir(cat(kernelopts(homedir),"/mapleprimes/phc")):

# Create an empty repository.
if FileTools:-Exists("phc.mla") then
  print("deleting old archive");
  FileTools:-Remove("phc.mla"):
end if;

       "deleting old archive"

LibraryTools:-Create("phc.mla"):
LibraryTools:-ShowContents( "phc.mla" );

                      []

read("phc.module"):
LibraryTools:-Save('phc', "phc.mla");
LibraryTools:-ShowContents( "phc.mla" );

    [["phc.m", [2019, 7, 2, 12, 59, 46], 41984, 1736]]

restart;
libname := cat(kernelopts(homedir),"/mapleprimes/phc"), libname:

with(phc);
    [cascade, computeResiduals, decompose, deflationStep,
     drawPaths, embed, eqnbyeqn, factor, filter,
     intersectWitnessSets, makeBertiniInput,
     makeBertiniStart, makeSolution, makeStartSystem,
     makeSystem, makeWitnessSet,
     printSolutions, printSystem, readSolutionsBertini,
     refine, removeBertiniFiles, setPHCloc,
     solutionsAppendToFile, solve, solveBertini,
     subsVariables, systemToFile, track, trackBertini]

setPHCloc(cat(kernelopts(homedir),"/mapleprimes/phc")):

T := makeSystem([x, y], [], [x^2+y^2-1, x^3+y^3-1]):

sols := solve(T):

printSolutions(T, sols);
(1) [x = -.67793e-31-.61630e-32*I, y = 1.0+.15407e-31*I]
(2) [x = -1.0+.70711*I, y = -1.0-.70711*I]
(3) [x = -.15407e-31-.11556e-31*I, y = 1.0+.13710e-31*I]
(4) [x = 1.0-.53357e-35*I, y = .79194e-30-.27271e-30*I]
(5) [x = -1.0-.70711*I, y = -1.0+.70711*I]
(6) [x = 1.0-.38519e-32*I, y = .30815e-32-.24652e-31*I]

The error message is coming from your attempts to invert Matrix V. But V is singular and noninvertible.

By the way you are attempting Matrix-Matrix multiplication, and not Vector-Matrix multiplication, since if the inverse of V were obtained then it would be a Matrix. V is a Matrix. There are no Vectors in your example, and VectorMatrixMultiply is not the right command.

You could also use the `.` command (just the dot), for multiplying Vectors and Matrices.

But your central problem is that V is noninvertible. There's a minor possibility that you could utilize a pseudo-inverse of V. I include it on the off chance. It may well not be what you want or expect.

Download matinv.mw