@mmcdara If you are using Maple 2015, or if the OP is using Maple 17, then it is significantly faster to do it like so:

Download Try_This_2_ac.mw 

In older versions there is (unnecessary, unwanted) wiping of cached information that occurs even when the parameters are set to their same, extant values. This causes considerable inefficiency, that can be avoid by the simple check that I show here. I believe that inefficiency is fixed in recent versions.

@ahmeng From my Maple 2017.2 for 64bit Linux,

Download isect.mw

When solve detects the explicit float coefficiant that I introduced into the example it tries to resolve the ensuing approximate RootOfs (via evalf, etc). So it's similar to evalf after allvalues. For me it is not fast on first attempt (and being fast upon repeats means nothing here -- it just remembers some results internally).

The maxsols for fsolve (univariate non-polynomial) usage can be traced through to a iterated use of Rootfinding:-Nextzero, which is yet another approach that one could do up manually.
   showstat(fsolve::scalarmultiple)
It is true that the maxsols option is not documented for the univariate nonpolynomial case. But it is documented for the univariate polynomial case. Your screenshot doesn't show even that, though, since it seems you are citing ?fsolve and not ?fsolve,details .

In your example, L is a table and not a list.

I changed this comment from another thread into its own Question. The subject matter is quite different.

Your Question lacks detail and clarity.

Your Matrix B has 22 rows and 4 columns.

What do you hope to get, when you write terms such as  B[1..1,[17, 18]] ? Do you understand that refers to the (nonexistent) 17th and 18th column entries of the first row?

What parts of B are you trying to plot?

How should we "correct" your code when you didn't tell us anything about what you wanted it to do?!

@dharr I adjusted your Answer slightly, and inlined images (right-click exporting the plots) so that it would be certain that Mapleprime's own rendering wasn't interfering.  I hope that's OK.

@nm You have G as  x*sin(Pi*x)  instead of  x*sin(Pi*y) .

But even if your second result were to attain, it has x restricted to integers and hence not free.

@dondomingo But you won't know home many solutions there are and how many names to use, in general. I don't think that is a good methodology.

I agree, nobody wants to have to re-type or copy&paste from output. It's not good programming to have to do so.

Just assign the full list of solutions to a single name, let's call it Sol, say.

Then you can get at the first solution by indexing, as Sol[1]. And then you evaluate any formula involving x and y using 2-argument eval (as showed in my Answer above). You could even access just the simplest formula x, or y.

Eg,
   eval( x, Sol[1] );
   eval( y, Sol[1] );
   nops(Sol);    # how many solutions there are
   eval(  x^2+sin(y-x),  Sol[3] );
etc.

@alijp I don't understand why you didn't show that when you initially asked the Question.

Now you've introduced the `i` parameter when you already had two (r and phi). So you don't need `i`, now.

Is your problem with the differentiation at symbolic order n?

(I have marked your Question as being for Maple 17.)

I could point out that a function might by undefined on an interval between discontinuities (including a pair returned by discont). And iscont could return true for the test on that open interval. This might matter to you.

@Anthrazit Your string collection omitted the opening and closing square brackets, needed so that after parsing a list would be passed to Join.

joinStr_ac.mw

This is an example of what I meant about confusing multiple concatenation, since all the little fiddly strings are wrapped in their own quotes. I think that an sprintf format is more legible and less error-prone.

@Thomas Richard Thanks for the detail.

I mentioned those two particular versions because I believe those are the two that the OP uses.

@Carl Love For the sake of completeness, I'll mention that the collapsed calling sequence for a multiple integral can be more than just a syntactic convenience (versus the iterated calling sequence).

When computing the value, the ranges of the outer variables can be used during the computation of the inner integrals, via assumptions on the relevant names of the outer variables of integration. It turns out that this can quite often be useful, and the collapsed form is used in some VectorCalculus and MultivariateCalculus in Maple 2020.

@mmcdara I think that I didn't understand your original goal. Based on the original code you gave (and with those R commands slightly in mind) I had the impression that you were trying to get a kernel-density-plot based on given Samples, rather than from symbolic RVs.

Hence the suggestions to try surfdata with dimension=2, and direct writing to an Image format, etc. I didn't mention plots:-listdensityplot, which has a possibly useful smooth=false choice but lacks the convenience of surfdata to specify the ranges for the axes. I was working on the presumption that one would be starting from (only) samples of unknown distributions. I suppose there may be room for both goals.

It might be interesting to tune up your scheme for using PDF and the symbolid RVs. As is, each call to KD2D calls the K procedure N times, and each K call generates its own procedure as return value, and each call to those returns call the f procedure twice, and each call to the f procedure generates its own procedure as return value. Alas I don't have any time for that fun these days.

I notice that you bumped down your Sample size to N=10^3, perhaps for time cost and perhaps also because at N=10^4 the points occlude the filled color gradient of the density plot. And they collide so much that the visual gradient is lost. (BTW, you could switch the order of plots in the display, so that at least the black lines of the contours lay in front of the possibly dense scatter-point rendering.)

So here is a quick little thing. In the case of a large Sample there is some visual benefit to plotting a scattering using points smaller than symbolsize=1. By replacing the points with very thin and short lines we can get a more visually distinguishable mist.

>  restart: >  with(plots): with(Statistics): >  f := (m, s) -> unapply(PDF(Normal(m, s), t), t): K := (mx, my, sx, sy) -> (x, y) -> f(mx, sx)(x) * f(my, sy)(y): >  N := 10^3: A := Sample(RandomVariable(BetaDistribution(2, 3)), N*60): B := Sample(RandomVariable(BetaDistribution(3, 2)), N*60): h := 0.07: KD2D := (x, y) -> add(K(A[n], B[n], h, h)(x, y), n=1..N)/N: >  CP := contourplot(KD2D(x, y), x=0..1, y=0..1, view=[0..1, 0..1],                   filled=true, coloring=["MediumBlue", "White"],                   contourlabels=false, thickness=0): >  SP := ScatterPlot(A, B, symbol=point, color="Gold", symbolsize=1): >  dat := op([1,1],SP): newSP := subsop([1,1]=Matrix(3*op([1,1],dat),op([1,2],dat),                              (i,j)->`if`(irem(i,3)=1,dat[(i+2)/3,j],                                          `if`(irem(i,3)=2,1e-3+dat[(i+1)/3,j],                                               Float(undefined)))),                 subsindets(SP,specfunc(anything,STYLE),                            ()->[STYLE(LINE),THICKNESS(2)][])): >  display(Array([display(SP, CP, size=[400,400]),                display(newSP, CP, size=[400,400])]));         >

Download 2D_density_plot_with_Kernel_ac0.mw

I made no changes for efficiency to the earlier code. And I would assemble the replacement Matrix of plotting data differently -- using BlockCopy rather than an initializer with nested `if` calls. And one has to be very careful with the scaling of the short, thin lines relative to the plot view, or else weird visual effects can occur. This is not production level code. I'm really just trying to show that I like the idea of symbolsize=0.1 and so on. In these situations the plot may look less like it was drawn with crayons.

I don't see any information in your attachment that could be used to restore anything.