Both answers clarify my doubts (thank you @acer), but @Carl Love 's one goes above and beyond with all the refreshers on the mathematical notions and the digressions. Hence the best answer.

@Carl Love thanks.

1. How to do this replacement within fsolve() command?
2. Since you brought up convergence...what's the standard (implicit) convergence of fsolve? How to output the error in case I want to have a better idea of the magnitudes compared to the numerical solutions given (apart from eval(equations,solutions))?

@tomleslie amazing answer! I am learning a lot. Thank you for all the details. I see that you provided fsolve() with specified starting points in all calibrations, not just those with "missing" values. Is this a good practice more in general? Should I always provide starting guesses or only if I encounter issues from some iteration onwards?

@Carl Love thank you for teaching me about this.

@Carl Love oh okay I understand now. Thanks!

@Carl Love the 1 plot per file solution (2nd option) seems like what I need. 

Since I have 8 consecutive plots:-display( seq( plot( [$ 20], series1of8[1..20,j], color=cols1[j]), j=1..3)); (one for each one of the 8 series), does it mean I first need to combine them all together as a plot sequence? (Is this what was implicit in Plots:= <seq(plot([$20], 'color'= cols1[j]), j= 1..3)>: ?)

@tomleslie in your first comment/answer you wrote "you have to ensure that all the data is available for all the plots whihc you want!"

I found missing data for data series 7 and 8: someplots_cal7_cal8_bounds_issue.mw

I found that series 7 stops at run 33, while series 8 at run 52. Of course if I only plot the first 20 elements or so there's no issue...but why do I have substantially less data for these two calibrations?

Separately, I noticed that your fsolve command for lam did not include the positive range for the lambdas, e.g. lam := fsolve(eval(MyEqs_cal, ncal3)) instead of lam := fsolve(eval(MyEqs_cal, ncal3), {lambda__1 = 0 .. infinity, lambda__2 = 0 .. infinity, lambda__3 = 0 .. infinity});...did you do this deliberately? For what reason?

@dharr thank you so much for the efforts. Interesting. I didn't think about it...I am learning a lot in this thread.

I was trying to look at a much simpler system: Analysis_quintic_solution.mw. Is there a way to find a closed-form real and positive solution to the quintic (with quadratic term missing) I find when I don't make any assumptions on the positivity of the two variables? Descartes rule of sign suggests that one such solution may be attainable, but how to find it explicitly? In the script are my attempts.

The system is unrelated to the one discussed over this thread, so feel free to ignore this last analysis, if you wish. Thank you for your help so far.

@dharr thank you for the exhaustive analysis. I checked a few times. EqN_rhou in the startup code seems to be the same I fed into the solver. However, as long as the solution solves the quadratic equation which was output by the solver, then it must also be a solution to the system...am I wrong?

I re-run that single calibration (took just 1 minute with the reduced 2-equations system). The quadratic equation produced (solver worksheet) is exactly the same as the one loaded in the analysis worksheet, but strangely the size of the solution differs slightly between the two worksheets. Is it the command sols_rhou := eval(SOLUTION_LAMBDAS_cal_rhou, '='tilde(%, {rho__u, lambda__2, lambda__3})) in the analysis worksheet that slightly reduces the length?

SOLVER WORKSHEET:

ANALYSIS WORKSHEET:

In any case, I reproduced the analysis calculations in the solver worksheets, directly on the solution output by the solver, and the lambdas still do not solve EqN_rhou. Why? How come that they solve the quadratic polynomial but not the original system?

@tomleslie thank you for the exhaustive answer !

@dharr understood. Do you mind double-checking whether I am missing anything from my analysis? 250523_SOLUTION_analytical.mw

I am not confident about it, I might have missed some corner solutions, as I got confused by the different domains...

From my (imperfect) analysis, it seems that for all 6 calibrations there's no solution allowing for both lambda_2 and lambda_3 to be positive. Would you confirm it? 

@tomleslie thank you very much for your clear answer. Yes, I realized the index issue but did not know how to collect each iteration in the right data structure. Your plots are very useful. I have two follow-up questions:

1. Can you show me calibration 12 and just one example representative for calibrations 9 to 11? These are correlations, so I'll need something like for INDEX from -1 by 0.002 to 1 do for calibrations 9 to 11 and for INDEX from -0.998 by 0.002 to 1 do for calibration 12 (since I noticed that from -1 gives me error for this calibration).
2. You wrote "You just have to decide exactly what you want to see". I wrote that in the body of my question. I'd like to use subplots. In particular, for the lambda plots I want to combine them for multiple calibrations (please refer to what I wrote above), while for the 12 beta plots I think is visually better to split the 6 data series so that I have just 2 for each sub-plot (but I do agree with you that there's no need for dualaxis plots).

@dharr thank you for the detailed analysis for gamma.

"if you look at RegularChains:-Triangularize that PolynomialSysyem is using, you can see how to more specifically add inequalities.."

• How about I just add the inequalities to the set of equations I pass to the solver (regardless of the solver I use)?

EDIT: I cannot do this apparently, since the only acceptable inputs are polynomials.

but it might be faster to process them after as needed."

• Does it mean that the quadratic solutions I obtain from the solver would be the same regardless of the assumptions on my variables and parameters (and only the roots to these quadratic equations depend on the assumptions)? 

@dharr a related question: is it possible that the absence of all-positive solutions is because the solver I used may ignore assumptions on the parameters (in the main script I defined them in the top execution block, same as I did in 220523_SOLUTION_reducedform_calibrations.mw)?

Because I know that standard solve() sometimes ignore these (sometimes even if you add the "assuming" option), but at least it outputs a warning...while I don't see such warning in the ouput of SolveTools:-PolynomialSystem. Does it mean that it took into account the assumptions while solving?

@dharr understood, thank you. I guess this is good news in the sense that from now onwards I can work with a 2-eqs-system. Do you mind sharing me your script for gamma? What about the other equations?