# How to transform a FindRoot in a NDSolve?

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:

tabexp=ParallelTable[{y1,y2/.FindRoot[Pinj-Fun[y1,y2],{y2,0.2,2},AccuracyGoal->25,PrecisionGoal->25,MaxIterations->2000]},{y1,0.2,10,0.9/200}];


This the NDSolve versione:

 sol=NDSolve[{Pinj==Fun[y1,y2[y1]],x'[y1]==1,x[0.2]==0.2,y2[0.2]==(y2y2/.FindRoot[Pinj-Fun[0.2,y2y2],{y2y2,2}])},y2,{y1,0.2,10}]


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;
Loss[y_,ye_]:=y5-ye5;
Pinj=0.01;


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.)