Plot the maximum of the real parts of the eigenvalues of a trancendental equation

Question

Basically I have this matrix

MatrixForm[{{-d (1/S), -S b, -S b, 0, 0, 0, 
   0, -S b, -S b}, {(d/2) (1/S - 1), -d - γ + b S, 0, 0, 0, 0, 0, 0, 
   S b}, {(d/2) (1/S - 1), 0, -d - γ + S b, 0, 0, 0, 0, S b, 
   0}, {0, γ (1 + G (lambda)), 0, -d, 0, 0, 0, 0, 0}, {0, 
   0, γ (1 + G (lambda)), 0, -d, 0, 0, 0, 
   0}, {0, -γ (G (lambda)), -(d + γ) + b S, 0, 
   0, -d - alpha phi ((d/(2 b)) (1/S - 1)), 0, -(d + γ) + b S, 
   0}, {0, -(d + γ) + b S, -γ (G (lambda)), 0, 0, 
   0, -d - alpha phi ((d/(2 b)) (1/S - 1)), 0, -(d + γ) + b S}, {0, 
   0, (d + γ) - b S, 0, 0, alpha phi ((d/(2 b)) (1/S - 1)), 0, -b S, 
   0}, {0, (d + γ) - b S, 0, 0, 0, 0, alpha phi ((d/(2 b)) (1/S - 1)), 
   0, -b S}}]

Every parameter has a value d=0.015, b = 180= alpha, Gamma = 52, S is a variable that depend on the paremters.

S is given by

S[phi_] := (-974001.6810000001` + 3.8452932405`*^6 phi - 
    Sqrt[9.486792745908259`*^11 - 9.247585725068096`*^11 phi + 
      8.244315219655941`*^12 phi^2])/(2 (-3.3705720000000005`*^6 + 
      3.3582924`*^6 phi))

And G(lambda) is a transcendental function

(E^(-1. λ)-11.910003349801029 +   E^(1. λ) (11.91 - 6. λ) - 
  5.910671637618375λ))/(0.015 + λ)^3

Thus for each value of phi in (0,4) , manually I calculate the value of S(phi) and then plot the eigenvalues like that

roots = Reduce[
  JacobisimetricoM[0.015, 52, 180, 1, 0.993, 180, 0.032, 
     0.27796141204714564, λ, (E^(-1. λ) (-11.910003349801029 + 
          E^(1. λ) (11.91 - 6. λ) - 
          5.910671637618375 λ))/(0.015 + λ)^3] == 0 && -1 < 
    Re[λ] < 70 && -3 < Im[λ] < 3, λ]

ListPlot[{Re[λ], Im[λ]} /. {ToRules[roots]}, 
 PlotLabel -> "Eigenvalues for ϕ=0.032 ", 
 PlotStyle -> {PointSize[0.02], Darker[Red]}, PlotRange -> All, 
 AxesOrigin -> {0, 0}, AxesLabel -> {"Real part", "Imaginary part"}]

where the equation JacobisimetricoM is the characteristic equation of the matrix (matriz form). (Lambda is the eigenvalue )

JacobisimetricoM[d_, γ_, b_, N_ , M_, alpha_, phi_, S_, λ_, Glambda_] := 1/(8 b^2 S^3) (d + λ)^2 (alpha d phi (N - S) (2 alpha d phi (N - S) (-d^2 + d (b S - 2 γ - λ) - γ (γ + \ λ) + b S (γ + Glambda γ + λ)) (d^3 N + λ (b S ((-1 + Glambda) γ - λ) + γ \ (γ + λ)) + d^2 (-b N S + S λ + N (2 γ + λ)) + d (N γ (γ + λ) + S λ (2 γ + λ) - b S (N (γ + λ) + S (-Glambda γ + λ)))) - (d + γ \+ λ) (alpha d phi (-N + S) - 2 b S (d + λ)) (b^2 d (N - S) S^2 (d + γ + Glambda γ + λ) - 2 (d + γ) λ (d + γ + λ) (d N +S λ) + b S (d^3 (-N + S) + 2 S λ^2 (γ + λ) + d^2 (N (-2 γ + λ) + S (2 γ + λ)) + d (S (γ^2 + γ λ - Glambda γ λ + 2 λ^2) + N (-γ^2 + γ λ + Glambda γ λ + 2 λ^2))))) + (d + γ + λ) \ (alpha d phi (-N + S) -2 b S (d + λ)) (b S (-d + b S - γ) (2 b S (d + λ) (b d (N - S) S + λ (d (N + S) + 2 S λ)) + alpha d phi (N - S) (d^2 (N + S) + 2 S λ (-b S + γ + λ) + d (-2 b S^2 + N (γ + λ) + S (γ + 3 λ)))) - (b S + λ) ((alpha d phi \(-N + S) - 2 b S (d + λ)) (b^2 d S^2 (-N + S) + 2 λ (d + γ + λ) (d N + S λ) - b S (d^2 (-N + S) + 2 S λ^2 + d (-N γ + S γ + 2 S λ))) - alpha d phi (N - S) (2 d^3 N + 2 S (b S - γ) (b S - γ - λ) \λ + d^2 (-b S (3 N + S) + 2 (2 N γ + N λ + S λ)) + d (b^2 S^2 (N + S) + 2 (N γ (γ + λ) + S λ (2 γ + λ)) - b S (N (3 γ + Glambda γ + 2 λ) +S (γ - Glambda γ + 4 λ)))))))

Although I obtain for each phi the eigenvalues, Now I would like to plot the maximum real part of eigenvalues as a function of phi. Anyone can help me?