some days ago I did a question about a problem with a findroot. The solution was to convert the FindRoot in a "fake" NDSolve. How can I do the same if the FindRoot has more equation ? I make here an example, this is the FindRoot version:


This the NDSolve versione:


This is more or less clear to me, the problem is that now i would like to solve in this way a FindRoot with two equations:

 tabexp=ParallelTable[{FindRoot[{Pinjx-Fun[y1,y2]- Loss[y1,1],Fun[y1,y2]-Loss[y2,1]},{y1,5},{y2,5},AccuracyGoal->25,PrecisionGoal->25,MaxIterations->2000]},{Pinjx,0.2,10,9.8/50}];

The rest of the code is here

 n[ω_?NumericQ,y_?NumericQ]:=(Exp[ ω/y]-1)-1;
Fun[y1_?NumericQ,y2_?NumericQ]:= NIntegrate[ω (n[ω,y1]-n[ω,y2]),{ω,0,10},MinRecursion->4,Method->{"GlobalAdaptive","MaxErrorIncreases"->100000,"SymbolicProcessing"->0,"SingularityHandler"->None},PrecisionGoal->5,WorkingPrecision->16];

1 Answer 1


Solving multiple equations as a DAE is similar to solving just one. Put them all in:

n[ω_?NumericQ, y_?NumericQ] := (Exp[ω/y] - 1) - 1;  (* sure this is right? *)
Fun[y1_?NumericQ, y2_?NumericQ] := 
  NIntegrate[ω (n[ω, y1] - n[ω, y2]), {ω, 
    0, 10}, Method -> {"GaussKronrodRule", "Points" -> 11}, 
   PrecisionGoal -> 10];
Loss[y_, ye_] := y - ye;
Pinj = 1/100;

Block[{Pinjx = 2/10},
 yy0 = {y1, y2} /.   (* initial starting values for y1, y2 *)
   FindRoot[{Pinjx - Fun[y1, y2] == Loss[y1, 1], 
     Fun[y1, y2] == Loss[y2, 1]}, {{y1, 2}, {y2, 3}}]

{ndsol} = NDSolve[
  {Pinjx - Fun[y1[Pinjx], y2[Pinjx]] == Loss[y1[Pinjx], 1], 
   Fun[y1[Pinjx], y2[Pinjx]] == Loss[y2[Pinjx], 1],
   {y1[2/10], y2[2/10]} == yy0,
   x'[Pinjx] == 1, x[2/10] == 2/10},
  {y1, y2}, {Pinjx, 2/10, 10}]

(I tweaked the NIntegrate code to make it converge more reliably & faster.)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.