@MaPal93 

I believe I do not have to answer your point 3 because the anser was in 141123_Problem_NoCorrelation_mmcdara_2.mw

About point 1

1. You can choose any value you want for lambda__2 and lambda__3 and thus get the corresponding value of lambda__1.
    
2. Is it true for all positive combinations of lambda_2 and lambda_3 that the roots of the numerators are also those of the denominators of my equations?
   I assume here you are talking about files 141123_Problem_NoCorrelation_mmcdara.mw and 141123_Problem_NoCorrelation_mmcdara_2.mw because I didn't evaluate the denominators in file 141123_Problem_mmcdara.mw.
   The answer to your underllined question is of course YES IT IS TRUE THAT ...
   Look to what I did (right below): where do you see that I used specific values for lambda__2 and lambda__3 that could give incidentlly nill denominators?

   sol_solve := solve(EqN, Vars):
   sol_solve := simplify(sol_solve, size):
   map(simplify, convert(eval(denom(Eqs_cal), sol_solve), list));
                              [0, 0, 0]
   

   sol_solve is a formal solution, not a numeric one.
   
   Look after the text The roots of the numerators are those of the denominators..., in this new file to see what happens when the lambdas are the roots of system EqN
   141123_Problem_NoCorrelation_mmcdara_3.
   
    

By the way, I've already spent hours to help you, don'y you think I deserve,at least, that you vote up?

@MaPal93 

Meaning of AllTotalDegrees:
AllTotalDegrees is a list of triples, if such a triple is [a, b, c]  then at least one of the equations to solve contains the term

​​​​​​(lambda__1)^a *  (lambda__2)^b *  (lambda__3)^c 

In the second script of yours (where you try to solve for the numerators):
It is you who tried to solve for the numerators, 
The lines

Eqs_redox := eval(Eqs_cal, p = .5); 
EqN := numer(Eqs_cal); nops(%); 
length(EqN); 
EqN_redox := eval(numer(Eqs_redox), lambda__1 = lambda__2); nops(%); 
length(EqN_redox)

are yours, not mine, I just added after them

sol_solve := solve(EqN, Vars):

Why did you use  SolveTools:-PolynomialSystem instead of solve?
(If I'm not mistaken it was a suggestion from @acer but I can't remember why he proposed that?)
Look to this new file from the text  Using SolveTools:-PolynomialSystem
141123_Problem_NoCorrelation_mmcdara_2.mw

Last advice: if you really want to understand what SolveTools:-PolynomialSystem is telling you begin to read these help pages

• RegularChains
• PolynomialRing
• ComprehensiveTriangularize
• PolynomialSystem

and read some papers or books about the subject... or act simpler.

@Carl Love 

Thanks for your feedback and, by the way, but don't feel obliged to spend too much time my code.

BTW, I wonder if I have not already ask a question about seaching all the occurrences of a specific entry in a hierarchical table, but  couldn't put the finger on it (maybe sevenor eight years ago)

@Carl Love 

Thanks your comments and corrections.

I tried writting a procedure that searches within a hierarchical table (a table whose entries can be tables, and so on) for all the indices where a given entry appears.
Here it is.
 

Download Recursive_search.mw

Explanation (highlighted is a wrong index, should be [6, 4])

here = [[4], [10]], [[0, -3]], [[0, a, r]], [[0, 6, 4]]

means that the target value (51) is an entry of T[4], T[0][-3], T[0][a][r], ...

Note that I used 

[entries(t,'pairs')]

instead of 

op(op(t))

because op(op(t)) is a list when procedure f is entered the first time, but a sequence when it is entered through thisproc.
(this seems weird to me)

If it's not too much trouble, could you help me correct this procedure, and suggest a more effective version?
Or do you think it would be better if I open a specific thread?

@Carl Love 

You're right, I didn't completely answered the question.
I edited my comment and completed it to account your remark.

@dharr 

My using isolate instead of solve is mainly a matter of laziness: isolate(equation, v) returns an equality that can directly be used in eval,  as solve (eq, v) requires writing v=solve (eq, v) to be used in eval.

@dharr 

I'm not at ease with eliminate and I often have to do several trials to get what I want.

Nice, I vote yup.

What is a(T) ?

@dharr 

Thanks

@shashi598

I adapted the Metropolis-Hastings' algorithm (developed to sample an unknown random variable) to find the values of alpha such that D(gamma)(xi__1) = 0.
(
More precisely, as the right boundary condition is gamma(xi__1)+xi__1* D(gamma)(xi__1) = 0, the previous condition is equivalent to gamma(xi__1) = 0, and it is equivalent to search for the values of alpha which make this identity true, or for values of alpha which minimize gamma(xi__1)^2.
I made this latter choice.
)

I got only 2 values (global minima of gamma(xi__1)^2) in the range alpha=1..10:

1. alpha = 3.7026096348364278322
2. alpha = 4.2985527068692070036

and 4 local minima (that is minima where gamma(xi__1)^2 is small but not null):

1. alpha = 5.3002873922249993805
2. alpha = 6.0544400684250000630
3. alpha = 7.5717744989250001680
4. alpha = 8.3303294592000002940

The type (local/global) of the minima are verified graphically.

You will note these values differ from those from the paper you cite.

Reply me if you still read this thread ad are interested to my code.

@Ronan 

Thank you Ronan

Could you check them?

@shashi598 

Explanation:

Is I'm not mistaken you want to find the value alpha such that diff(gamma(xi), xi) = 0 at xi=xi__1.
Which, give, that gamma(xi) + xi * diff(gamma(xi), xi) = 0 at xi=xi__1, is equivalent to find the value alpha such that gamma(xi) = 0 at xi=xi__1.

Here is the plot, with a logarithmic vertical axis, of gamma(xi__1)^2 in the range 1..6:


Then a focus on the range 3.5..4


and another on the range 4..6

Here are the values of alpha that minimize H(xi) = amma(xi__1)^2 when alpha is assumed to move in diffferent ranges:

This last plot completes the previous one by adding, in gray", the basins of attraction of red points, and in green those of the blue points.
Depending on your starting point in the range 4..6, NLPSolve will get anyone of the 5 points.
By default the starting point (if not specified with the option initialpoint) is the mid point of the searching range, thus 5 here.
Starting from 5 will return the first leftmost red point on this figure (that is your alpha2_result point).

Using DirectSearch may enable jumping out of this basin of attraction (the leftmost green rectangle) for it is not very deep.
Unfortunately the algorithm ends hitting the right boundary of the search domain (where the value of gamma(xi__1)^2 is lower than those at the red point, so it indeed minimizes the objective, at least in some way), but it misses the blue points because of their thinness and a bad tunning (I guess so) of the algorithm. 

Isn't it this a duplicate of thar previous question?

BTW, why do you persist in doing the same mistake (unknown=f(x) and parameter mu[f])?
Did you read my answer?


More of this your 4 equations are both identical except for the name of the dependent function, so their solutions will be fundamentaly identical it it is stupid to dsolve 4 equations.
The generic form of these equations is given by coeff(HPMeq, p, 1).
Before doing complex things  ask yourself if this equation can be solved formally; 

# example
test := coeff(HPMEq, p, 0);

# A solution is returned here
Without_bc := dsolve(test);   

# But not here
With_bc := dsolve({test, F[0](1) = 1, F[0](-1) = -1, (D(F[0]))(1) = 0, (D(F[0]))(-1)})

So trying to take the rhs of this last ouput generates the error message you got.

The only way is to dsolve numerically

{test, F[0](1) = 1, F[0](-1) = -1, (D(F[0]))(1) = 0, (D(F[0]))(-1)}

which implies to replace the many parameters by their numerical values.

Make an effort to understand all this and to submit us something that is worth us spending time answering, otherwise it's all a waste of  time for the people here.

@Saha 

 

Download eq1_only.mw