@acer thanks for describing the inner workings of fsolve().

You write: "Moreover (and I know this sounds facile, but I'm being serious), it can help to pause and try to define the notion of a "root" -- in this float context. Once you have such, you can program for it."
I don't understand what you have in mind.

Also "Using black-box procedures can separate the working precision for the functional evaluations (Digits now set withing the scope of the procedures) from the accuracy requirement in the residual check." So how to use operator-form instead of expressions, for the equations? What exactly does it mean? Are you confident that it's the residual check at INDEX~0.554 the root cause of the singularity?

@mmcdara thanks. 

You write: "Maybe reducing the step of INDEX before reaching this value or increasing the number of Digits would prove I'm wrong?"
How to reduce the step of INDEX dynamically? Are you thinking of a for sub-loop within the main for loop from 0.8731 onwards?

Moreover, I set Digits:=30 inside the "if decrease then", right before lam_ranges:=...
Since I am not seeing any improvement with that, I assume that placing the Digits command there is not what you had in mind...

Anyway, I understand if you are done looking at this thread and I thank you for your time and ideas so far. Actually it's already an improvement compared to stopping at the 600-something iteration and is nice to see that the discontintuity is actually not there. Perhaps I can think some more about it. Initially I thought that decreasing the lambda__i1 lower bound even more could be helpful but I think I now understand that it could simply be that there are no positive lambda__i1 (however small they are) beyond that threshold iteration...Maybe the dynamic adjustement of the step for INDEX is the technical solution we are looking for... Let's see if anyone else has other ideas.

UPDATE: See discontinuity_persists.mw (leaving aside what happens for INDEX>0.5):

1. Did your script contain a spelling mistake? See the two edited lines of code highlighted in electric blue
2. When plotting for INDEX from 1e-3 to 0.5, why the discontinuity at ~0.27 DOES NOT go away despite reducing the step to 1e-4? (despite plotting for INDEX from 0.250 by 1e-4 to 0.300 makes the discontinuity disappear...I thought that once the step is reduced to 1e-4 there wouldn't be any discontinuity anymore regardless of the plotting range)

@Carl Love thank you for entering the thread. You can take a look at this worksheet: fsolve_loop_breaks.mw.

Hopefully your fix tackles both problematic iterations, that is, the one breaking fsolve() as well as the one causing the discontinuity. It would also be useful to include some basic statistics about the error magnitudes over the whole loop. 

@mmcdara lambda_d1 = lambda_d2 and  lambda_i1 = lambda_i2 only if sigma__v1=sigma__v2. All 6 lambda are distinguished if these two parameters are set different (see for example fsolve_help_mmcdara_ncal1.mw).
By the way, what's the purpose of wrapping the array construction into the if lam::set then statement? What if you didn't do that? 

Two issues related to running your script for ncal2: fsolve_help_mmcdara_ncal2_breaks.mw

1. A discontinuity emerges for sigma__delta1~0.27
2. fsolve() "breaks" at the 614th iteration

How to avoid/prevent the discontinuity and how to "skip" the 614th iteration so that I can "fsolve() all the way" to the 1000th and final iteration and fully fill up the lamRes1 array?

@Axel Vogt thank you but as I replied to @mmcdara I don't think the negativity is the reason why fsolve() breaks after a specific iteration. The fact that fsolve is inside the loop is an important aspect of the issue I am trying to solve.

@tomleslie once wrote:
"When data 'goes missing' is because the fsolve() command is failing to find  solution for certain loop index values, at which point all subsequent code 'breaks'. In order to assist  the fsolve() command as much as possible, there are a couple of possibilities

1. one strategy is to restrict the region over which it has to search for a solution
2. second strategy is to provide fsolve with a "starting point" which is "close" to the desired answer"

My worksheet fsolve_does_not_work.mw uses (2) above, where the starting point for the i-th loop iteration is the result of the (i-1)-th loop iteration. The initial guess for the first iteration in each calibration, I took as {lambda__11 = 0.25, lambda__12 = 0.25, lambda__13 = 0.25, lambda__21 = 0.25, lambda__22 = 0.25, lambda__23 = 0.25, lambda__31 = 0.25, lambda__32 = 0.25, lambda__33 = 0.25}. This strategy generates values only up to iteration number 55 out of 1000.

@Carl Love suggested to replace fsolve non-answers with undefined:
"Like many things in numerical analysis, fsolve is an iterative process attempting to create a numeric sequence (or a sequence of numeric vectors) that converges in some sense. Now, I've learned over the years that despite our best efforts (e.g., using all the tricks that Tom has provided), there will always be some cases where fsolve doesn't converge. Indeed, I don't think that I've ever had a practical system of nonlinear equations with parameters and more than 50 sets of parameter instantiations where I could get fsolve to give an answer for all instantiations. So, whenever you get a non-answer from fsolve, whatever values that would've gone into your plots had it returned numeric values should be replaced by the keyword undefined."

My questions then are:

1. How to do this replacement within fsolve() command? Could I then replace such "undefined" with a simple interpolation of the last valid iteration and the next valid iteration? By doing so perhaps I could nicely plot all my curves up to the 1000th iteration.
2. What's the standard (implicit) convergence of fsolve? How to output the error within fsolve() command in case I want to have a better idea of the error magnitudes compared to the numerical solutions at each iteration?

@acer mmm, okay thank you. No other way to simply show a longer line symbol in the legend box? Even though I really didn't want to, your example actually made me realize that perhaps using three different colors would look clearer overall...

@mmcdara thank you for showing that perhaps my example was a bad one since you actually found that there are no solutions assuming positivity of all 6 lambdas (I think this is consistent with the fact that fsolve() loop breaks at the first iteration).

My question was different though: I would like to address fsolve() "fragility", i.e., I want to know how to help fsolve() to pin down solutions at each iteration.

Allow me to consider a different example. Here fsolve_works.mw fsolve() correctly generates the full array for 1000 iterations. Note that sigma_v1=sigma_v2=0.1. However, here fsolve_does_not_work.mw fsolve() correctly generates the full array for the first 55 iterations and 'breaks' after that. Note that everything is the same as the working worksheet except that sigma_v1=0.1 but sigma_v2=0.11 (they are only slightly different now). The fact that the first 55 iterations go through suggests that there are strictly positive solutions for all 9 lambdas (unlike my badly chosen original example), and yet fsolve() 'breaks' at a specific iteration. I want to know how to address this (maybe skip such problematic iterations and interpolate to fill these 'gaps'? How to do that? What else to do?)

@mmcdara @acer ? Would you please follow up? I would prefer to have my remaining doubts clarified in this thread instead of opening a new question...

@dharr thank you!! That's perfect. Yes the scaling is super easy to see now. Most likely will be this easy also for all the other functions. Moreover, the "new" notation and related reduction of numerical evaluations makes the deriv much more compact.

@dharr thanks a lot. Works perfectly for all lambda_1 and lambda_2 derivatives.

However, as you anticipated, it's harder to manually scale derivatives of functions of lambda_1 and lambda_2. See for example d(beta_1)/d(sigma_d). An old answer of yours could perhaps help you think of a procedure to automatically do such scaling. See the weblink I attached at the bottom of this worksheet:

(How come there is no built-in Maple command that returns something similar to indets() but with powers?)

EDIT: after some manual handling I found that the transformation deriv/(sigma__d*gamma) nicely scales this specific derivative, so it's all good for this one. What I did: 1) first evaluate deriv for sigma__v1=Gamma__1/sigma__d/gamma,sigma__v2=Gamma__2/sigma__d/gamma), 2) then focus on the first term only of deriv. I am not sure if I can do the same for the other derivatives and with similar ease...

Download deriv_autoscaling.mw

@acer thank you. I didn't know about RootFinding:-Parametric.

1. Why the plot doesn't change if I also include -1<rho and rho<1 among the equations to solve? I understood that regions 1, 4 and 5 have 0 solutions anyway but if I add the constraints on rho I should expect 2 regions instead of 5, right?
2. I don't understand SampleSolutions(m,2)=0.56 and SampleSolutions(m,3)=0.51 (even after reading help page). How are these numeric values found?
3. Looking at the big picture, how to reconcile this CellPlot with A) my plot3d of Eq and B) @mmcdara Sturm's analysis. I am having a hard time putting all together.

Thanks again for following up on this thread.

@mmcdara I am finally back to this. Since it's my fault not having specified Gamma>0 as a pre-requisite of this question, I "fixed" your analysis worksheet myself to consider only Gamma>0: rho-analysis_mmcdara_Gammapositive_MaPal.mw.

1. Would you please double check whether I did everything correctly? I highlighted in light blue my changes.
2. What's the conclusion of Sturm's analysis? My understanding: for rho=0.6 there are 4 distinct real roots when 0<Gamma<max(4 roots-green bar), 2 real roots when Gamma is in yellow bar, and 2 real roots when Gamma is in red bar. Am I correct?
3. We now confirmed the existence of 4 real Lambda roots when Gamma values are in the green bar. How to isolate them? That is, what is the closed form of these 4 roots? How do I formally prove that just one out of four is positive (since this is what plot3d(Eq) for Gamma>0 and -1<rho<+1 suggests, see plot attached below)? 
4. What about the two real roots in the yellow rectangle and the two real roots in the red rectangle? Why plot3d shows just 4 roots for any Gamma in the range 0<Gamma<10? 
   Thank you for addressing my remaining doubts @mmcdara .

@acer sure, no problem. Or perhaps there's an easier way to infer the signs of these partial derivatives given that I know the signs of l1, l2 in their specified form and of their partial derivatives wrt to their variables Gamma_1 and Gamma_2 (thanks to @dharr's nondimGamma1Gamma2.mw, see the plots at the bottom). 

Note also that Gamma_1=gamma*sigma_d*sigma_v1>0 and Gamma_2=gamma*sigma_d*sigma_v2>0 because all the underlying parameters are strictly positive. These are the partial derivatives for L1, L2, B1, B2, A1, A2, A2S and A1S when l1 and l2 are left unspecified:

Download derivatives_signs.mw

@dharr and @acer thank you for your comments. Do you agree that both approaches are equivalent?

Moreover, doubling the digits doesn't solve the jumps...I also just noticed that the partial derivatives of L1 and L2, while still convoluted, are significantly simpler than those for B1, B2, A1, A2, A2S and A1S. Maple takes forever to find the derivatives for the latter 6 functions (the diff() step right before plot3d(preprocess())...). For example, figuring out just dB1dSd and its 3D plot is taking Maple more than 45 minutes...

So perhaps "manipulation to a nicer form" for both l1 (the RootOf()) and l2 (rational function of greatest l1 root) is needed here, but I understand that it is not easy. @dharr you mention continued fractions...is it virtually equivalent to an approximation that is accurate for an upper bounded range of positive Gamma_1 and Gamma_2 values? @acer do you have some ideas on how to manipulate or, if easier, directly approximate the bivariates l1 and l2 (and measure the accuracy)?  

If needed, I can spin off the question about the manipulation/approximation of l1 and l2 into a separate question (since it's not about non-dimensionalization anymore) but moderators probably don't want that.

@dharr that makes sense, it's just a scaling. Just replacing T with Gamma__2/Gamma__2 confirmed what you say. Thank you.

In a comment above you wrote:
"The problem may be close roots not being properly resolved by the default numerical routines used by plot3d. Also, the order of the RootOfs if there is more than one positive root is to have the smallest one first. The solution is to be careful about always choosing the most positive root for Lambda__1."

I think a similar issue emerges when I try to plot3d partial derivatives of some functions of Lambda__1(Gamma__1,Gamma__2)=RootOf(degree-10-polynomial) and Lambda__2(Gamma__1,Gamma__2) (where Lambda__2 should be interpreted as a rational function of the greatest root of Lambda__1 as in your nondimGamma1Gamma2.mw above):
derivatives_jumps.mw

In the attached script derivatives_jumps.mw you see jumps/discontinuities emerging in the derivatives of L1 and L2, seemingly for Gamma__1=Gamma__2~1...I am not sure whether these jumps really exist though, given that all parameters are strictly positive and in nondimGamma1Gamma2.mw above you found the partial derivatives of Lambda__1 and Lambda__2 wrt to Gamma__1 and Gamma__2 to be either strictly positive of strictly negative...

then how to have a "better"/smoother numerical evaluation in these more complicated partial derivative cases and for the ranges Gamma__1=0.01..5,Gamma__2=0.01..5? Quoting you, I know I should "make a function that reliably returns the numeric value of the largest root of Lambda_1 first and then change the RootOf form just to Lambda__1 (similar to if we had solved without backsubstitution) and use the numerical solution for Lambda__1 to find Lambda__2" but where to exactly implement this in derivatives_jumps.mw without messing up the function assignments? Perhaps when I specify the forms for l1 and l2 and their respective aliases? Moreover, for the derivatives I guess I would need the denoms of Lambda__1 and Lambda__2...