PE:=PowerExpand
FS:=FullSimplify

(* do 1-loop and 2-loop ourselves *)
SetAttributes[δ,Orderless]
Clear[i,j,k,l,m,n,o,p,a]
Unprotect[Times,Power];
δ[i_,j_]δ[j_,l_]:=δ[i,l]
ϕ[i_] δ[i_,j_]:=ϕ[j]
((ϕ[j_])^2):=ϕ2
δ[k_,k_]:=n
((δ[k_,l_])^2):=n
Protect[Times,Power];
int[i_,j_,k_,l_]:=(δ[i,j]δ[k,l]+δ[i,k]δ[j,l]+δ[i,l]δ[k,j] )/3
Print["1-loop tadpole"]
ϕ[i]ϕ[j]int[i,j,k,k]//Expand//FS
Print["1-loop banana"]
ϕ[i]ϕ[j]int[i,j,k,l]int[k,l,o,p]ϕ[o]ϕ[p]//Expand//FS
comb1=1/4! (3!)^2 2/2!
Print["total 1-loop contribution:"]
banana=%% %%%
Print["2-loop banana"]
ϕ[i]ϕ[j]int[i,j,k,l]int[k,l,o,p]int[o,p,q,r]ϕ[q]ϕ[r]//Expand//FS
combbanana2=(1/4!)^2 *3 (3!)^3* 2^2/3!
banana2=%% %//FS;
Print["2-loop LH"]
ϕ[i]ϕ[p]ϕ[q] ϕ[r]int[i,j,k,l]int[k,l,o,p]int[j,o,q,r]//Expand//FS
comb2LH=2(1/4!)^2 *(3!)(3 4)^2 
LH=%% %;
Print["2-loop WF"]
ϕ[i]ϕ[j]int[i,j,k,k]//Expand//FS
combWF2=(4/4!)^2  3!
WF2=%% %;
Print["sums"]
banana2/ϵ^2 /((ϕ2)^2)+ LH/((ϕ2)^2)(1/2/ϵ^2+1/4/ϵ)-(banana/((ϕ2)^2)/ϵ)^2-2 WF2/ϕ2 ( 1/(8 ϵ))
Series[%,{ϵ,0,1}]//FS
%/.a->1//Normal//FS
Print["seem to have different norms, multiply by 4, and by 2 for 2-loop contribution"]
%% 4 2

βuou=((8 u)/3)-((14 (u^2))/3)+2 v-((22 u v)/3)-((5 (v^2))/3)+7 (v^3)-ϵ+u (v^2) ((107/4)+8 ζ[3])+(u^3) ((370/27)+((88 ζ[3])/9))+(u^2) v ((659/18)+((64 ζ[3])/3))+(v^4) (-(477/16)-3 ζ[3]-6 ζ[4])+(u^3) v (-(15967/81)-((4856 ζ[3])/27)+((340 ζ[4])/9)-((2560 ζ[5])/9))+(u^2) (v^2) (-(13433/54)-((1456 ζ[3])/9)+((64 ζ[4])/3)-((2000 ζ[5])/9))+(u^4) (-(24581/486)-((4664 ζ[3])/81)+((352 ζ[4])/27)-((2480 ζ[5])/27))+u (v^3) (-(4867/36)-50 ζ[3]-8 ζ[4]-((160 ζ[5])/3))+(v^5) ((158849/1152)+((1519 ζ[3])/24)-18 ((ζ[3])^2)+((65 ζ[4])/2)+2 ζ[5]+75 ζ[6])+u (v^4) ((1336801/1728)+((5495 ζ[3])/12)-((190 ((ζ[3])^2))/3)+((141 ζ[4])/2)+((1145 ζ[5])/3)+((575 ζ[6])/3)+441 ζ[7])+(u^5) ((17158/81)+((27382 ζ[3])/81)+((1088 ((ζ[3])^2))/27)-((880 ζ[4])/9)+((55028 ζ[5])/81)-((6200 ζ[6])/27)+((25774 ζ[7])/27))+(u^2) (v^3) ((2281727/1296)+((37789 ζ[3])/27)-((544 ((ζ[3])^2))/9)-((337 ζ[4])/3)+((17444 ζ[5])/9)-((1600 ζ[6])/9)+2352 ζ[7])+(u^4) v ((537437/486)+((116759 ζ[3])/81)+((3148 ((ζ[3])^2))/27)-((10177 ζ[4])/27)+((75236 ζ[5])/27)-((24050 ζ[6])/27)+((11564 ζ[7])/3))+(u^3) (v^2) ((1314497/648)+((171533 ζ[3])/81)+((1384 ((ζ[3])^2))/27)-((23105 ζ[4])/54)+((96794 ζ[5])/27)-((25400 ζ[6])/27)+((14210 ζ[7])/3));

βvov=4 u-((82 (u^2))/9)+3 v-((46 u v)/3)-((17 (v^2))/3)-ϵ+(v^3) ((145/8)+12 ζ[3])+(u^3) ((821/27)+((224 ζ[3])/9))+u (v^2) ((131/2)+48 ζ[3])+(u^2) v ((325/4)+64 ζ[3])+(u^2) (v^2) (-(10661/18)-724 ζ[3]+184 ζ[4]-((3440 ζ[5])/3))+(u^3) v (-(12349/27)-((5312 ζ[3])/9)+((440 ζ[4])/3)-960 ζ[5])+u (v^3) (-(1004/3)-387 ζ[3]+96 ζ[4]-600 ζ[5])+(u^4) (-(19679/162)-168 ζ[3]+40 ζ[4]-((7280 ζ[5])/27))+(v^4) (-(3499/48)-78 ζ[3]+18 ζ[4]-120 ζ[5])+(v^5) ((764621/2304)+((7965 ζ[3])/16)+45 ((ζ[3])^2)-((1189 ζ[4])/8)+987 ζ[5]-((675 ζ[6])/2)+1323 ζ[7])+(u^5) ((389095/729)+((259358 ζ[3])/243)+((13288 ((ζ[3])^2))/81)-((27166 ζ[4])/81)+((179696 ζ[5])/81)-((63500 ζ[6])/81)+((28420 ζ[7])/9))+u (v^4) ((1067507/576)+((35083 ζ[3])/12)+288 ((ζ[3])^2)-((3697 ζ[4])/4)+5920 ζ[5]-2100 ζ[6]+7938 ζ[7])+(u^4) v ((1279979/486)+((784621 ζ[3])/162)+((18154 ((ζ[3])^2))/27)-((83837 ζ[4])/54)+((275510 ζ[5])/27)-((98975 ζ[6])/27)+((43120 ζ[7])/3))+(u^2) (v^3) ((3633377/864)+((125459 ζ[3])/18)+((2266 ((ζ[3])^2))/3)-2263 ζ[4]+14328 ζ[5]-((15575 ζ[6])/3)+19404 ζ[7])+(u^3) (v^2) ((9309907/1944)+((224804 ζ[3])/27)+((3032 ((ζ[3])^2))/3)-((73018 ζ[4])/27)+((155692 ζ[5])/9)-6300 ζ[6]+23912 ζ[7]);

η=((u^2)/9)-((2 (u^3))/27)+((125 (u^4))/324)+((u v)/3)-(((u^2) v)/3)+((125 (u^3) v)/54)+((v^2)/6)-((3 u (v^2))/8)+((145 (u^2) (v^2))/36)-((v^3)/8)+((65 u (v^3))/24)+((65 (v^4))/96)+(u^2) (v^3) (-(13741/432)+((67 ζ[3])/18)-10 ζ[4])+(u^3) (v^2) (-(58177/1944)+((191 ζ[3])/54)-((260 ζ[4])/27))+u (v^4) (-(18545/1152)+((15 ζ[3])/8)-5 ζ[4])+(u^4) v (-(3010/243)+((115 ζ[3])/81)-((110 ζ[4])/27))+(v^5) (-(3709/1152)+((3 ζ[3])/8)-ζ[4])+(u^5) (-(1204/729)+((46 ζ[3])/243)-((44 ζ[4])/81));

νinverse=2-((2 u)/3)+((5 (u^2))/9)-((37 (u^3))/18)-v+((5 u v)/3)-((37 (u^2) v)/4)+((5 (v^2))/6)-((251 u (v^2))/24)-((7 (v^3))/2)-(u^2) (v^2) (-(9199/108)-((26 ζ[3])/3)-((52 ζ[4])/3))-u (v^3) (-(4243/72)-((19 ζ[3])/3)-12 ζ[4])-(u^3) v (-(7765/162)-((136 ζ[3])/27)-((88 ζ[4])/9))-(v^4) (-(477/32)-((3 ζ[3])/2)-3 ζ[4])-(u^4) (-(7765/972)-((68 ζ[3])/81)-((44 ζ[4])/27));

order=2;
βuou/.{u-> κ u,v->κ v};
Series[%,{κ,0,order}];
Normal[%]/.κ->1;
Coefficient[%,u^2 v^2]/.{u->0,v->0}//Expand;
(* checked up to 4-loop order with paper *)

βvov/.{u-> κ u,v->κ v};
Series[%,{κ,0,order}];
Normal[%]/.κ->1;
Coefficient[%,u^3 v]/.{u->0,v->0}//Expand;
(* checked up to 4-loop order with paper *)

η/.{u-> κ u,v->κ v};
Series[%,{κ,0,order}];
Normal[%]/.κ->1;
Coefficient[%,u^3 v]/.{u->0,v->0}//Expand;
(* checked up to 4-loop order with paper *)

νinverse/.{u-> κ u,v->κ v};
Series[%,{κ,0,order}];
Normal[%]/.κ->1;
Coefficient[%,u^0 v^4]/.{u->0,v->0}//Expand;
(* checked *)

fixedpoint={ϵ->δ^2,u-> Sum[δ^n a[n],{n,1,order}],v-> Sum[δ^n b[n],{n,1,order}]};

ββ={βuou/.fixedpoint,βvov/.fixedpoint};
gen[n_]:={SeriesCoefficient[ββ[[1]],{δ,0,n}]==0,SeriesCoefficient[ββ[[2]],{δ,0,n}]==0}

(* find the δ = sqrt[ϵ] expansion for the two solutions, whichsol=1,2 *)
For[whichsol=1, whichsol≤2,whichsol++,
sol={};
sol1=gen[1]//Solve//Flatten;
X=gen[2]/.sol1//FS;
sol2=Solve[X];
sol={sol1,sol2[[whichsol]]}//Flatten;
X=gen[3]//.sol//FS;
sol3=Solve[X,{a[2],b[3]}]//Flatten;
sol={sol,sol3}//Flatten;
X=gen[4]//.sol//Simplify;
sol4=Solve[X,{a[3],b[4]}]//Flatten;
sol={sol,sol4}//Flatten;
X=gen[5]//.sol//Expand//Simplify;
sol5=Solve[X,{a[4],b[5]}]//Flatten;
sol={sol,sol5}//Flatten;
solution[whichsol]=sol;
]

(* check that the found solution is a fixed point *)
{βuou,βvov}/.fixedpoint//.sol;
Series[%,{δ,0,order}]
{βuou,βvov}/.fixedpoint//.sol;
Series[%,{δ,0,order}]

(* give the fixed point *)
{u,v}/.fixedpoint//.sol;
Series[%,{δ,0,order-1}]

(* the two solutions give different exponents *)
η/.fixedpoint//.solution[1];
η1=Series[%,{δ,0,order-1}]//FS
η/.fixedpoint//.solution[2];
η2=Series[%,{δ,0,order-1}]//FS
η1-η2

num:={ζ[n_]:> N[Zeta[n]]};

(* stability of the Gaussian fixed point *)
{{D[βuou u ,u],D[βuou u,v]},{D[βvov v ,u],D[βvov v,v]}}//.{u->0,v->0};
Series[%,{δ,0,2}]//Normal//MatrixForm
Eigenvalues[%];
Series[%,{δ,0,2}]/.num

(* this is the fixed point for v=0 *)
Series[βuou+ϵ/.v->0,{u,0,order}];
fp1u=InverseSeries[%]/.u->ϵ;
(* check *)
βuou/.{v->0,u->fp1u}
(* stability of this fixed point *)
{{D[βuou u ,u],D[βuou u,v]},{D[βvov v ,u],D[βvov v,v]}}//.{u->fp1u,v->0};
Eigenvalues[%];
Series[%,{ϵ,0,order}]/.num
(* so this one is completely stable *)

(* this is the fixed point for u=0 *)
Series[βvov+ϵ/.u->0,{v,0,order}];
fp2v=InverseSeries[%]/.v->ϵ;
(* check *)
βvov/.{v->fp2v,u->0}
(* stability of this fixed point *)
{{D[βuou u ,u],D[βuou u,v]},{D[βvov v ,u],D[βvov v,v]}}//.{u->0,v->fp2v};
Eigenvalues[%];
Series[%,{ϵ,0,order}]/.num
(* so this one has one stable and one unstable direction *)

(* stability of the non-trivial fixed point with u,v =!= 0 *)
{{D[βuou u ,u],D[βuou u,v]},{D[βvov v ,u],D[βvov v,v]}}/.fixedpoint//.solution[1];
Series[%,{δ,0,order-1}]//Normal;
Eigenvalues[%];
Series[%,{δ,0,order-1}]/.num

(* plot all the fixed points *)
epsilon=0.25;
order=5;
(* choose the region you want to see *)
region=1;
Which[region==1,vrange={-0.15,0.55};
urange={-0.4,0.25},
region==2,
vrange={-1,0.5};urange={-0.4,1},
region==3,vrange={-0.15,2};urange={-0.5,2}];
δv=vrange[[2]]-vrange[[1]];
δu=urange[[2]]-urange[[1]];
Print["u in ",urange]
Print["v in ",vrange]
(* get the series exansions of all four fixed points; they are plotted in blue *)
{v,u}/.fixedpoint//.solution[1]/.a[_]->0/.num;
Series[%,{δ,0,order-1}]//Normal;
%/.δ->Sqrt[epsilon];
s1=Series[{{0,0},{0,fp1u},{fp2v,0}},{ϵ,0,order}]/.num//Normal;
{v,u}/.fixedpoint//.solution[2]/.a[_]->0/.num;
Series[%,{δ,0,order-1}]//Normal;
%/.δ->Sqrt[epsilon];
s2=Series[{{0,0},{0,fp1u},{fp2v,0}},{ϵ,0,order}]/.num//Normal;
Join[s1,s2,{%%}]/.ϵ->epsilon;
LPblue=ListPlot[%,PlotStyle->Blue];
(* sample the region 100 times to find all fixed points, they will be plotted in red *)
{βvov v,βuou u}/.num/.ϵ->epsilon;
ser=(Series[%/.{u->κ u,v->κ v},{κ,0,order+1}]//Normal)/.κ->1;
Off[FindRoot::lstol];Off[FindRoot::cvmit];
Table[{v,u}/.FindRoot[ser,{v,vrange[[1]]+δv  RandomReal[]},{u,urange[[1]]+RandomReal[]δu},MaxIterations->1000],{100}]//Chop;
On[FindRoot::lstol];On[FindRoot::cvmit];
Round[%%,10^(-6)]//Union;  
(* rounded, so that fixed points are not multiply found *)
LP=ListPlot[%,PlotStyle->Red];
(* starting points for the flow *)
pts=Table[{vrange[[1]]+RandomReal[]δv,urange[[1]]+RandomReal[]δu},{300}];
{-βvov v,-βuou u}/.num/.ϵ->epsilon;
sf=(Series[%/.{u->κ u,v->κ v},{κ,0,order+1}]//Normal)/.κ->1;
streamColor[x_]=Cyan;
SP1=StreamPlot[sf,{v,vrange[[1]],vrange[[2]]},{u,urange[[1]],urange[[2]]},PerformanceGoal->"Quality",AspectRatio->Automatic,StreamPoints->{pts,Automatic,Scaled[2]},StreamScale->None,StreamColorFunction->streamColor];
SP2=StreamPlot[sf,{v,vrange[[1]],vrange[[2]]},{u,urange[[1]],urange[[2]]}];
Show[LPblue,SP1,SP2,LPblue,LP,PlotRange->{vrange,urange},AspectRatio->Automatic,AxesLabel->{v,u}]