Inv = i0 - i1*r[t];
S = (1 - b)*(y[t] - T) - a;
IS = (Inv + G == S + T);
ISnerovnovaha = Inv + G - (S + T);
L = C0 + C1*y[t] - C2*r[t];
LM = (L == M);
LMnerovnovaha = (L - M);
y'[t] == 
      Collect[Expand[k1*ISnerovnovaha], {y[t], 
          r[t]}]
r'[t] == 
      Collect[Expand[k2*LMnerovnovaha], {y[t], r[t]}]

C0 = 356.6;
C1 = 0.4;
C2 = 18;
M = 357.1;
i0 = 150.66;
i1 = 10.24;
a = 5;
b = 0.8826;
T = 136.51;
G = 171.43;
k1 = 7;
k2 = 0.5;
y'[t] == 
      Collect[Expand[k1*ISnerovnovaha], {y[t], 
          r[t]}]
r'[t] == 
      Collect[Expand[k2*LMnerovnovaha], {y[t], r[t]}]

Solve[{LM, IS}, {y[t], r[t]}]

A = {{Coefficient[Simplify[k1*ISnerovnovaha], y[t]], 
            Coefficient[Simplify[k1*ISnerovnovaha], r[t]]}, {Coefficient[
              Simplify[k2*LMnerovnovaha], y[t]], 
            Coefficient[Simplify[k2*LMnerovnovaha], r[t]]}};

MatrixForm[A]

Eigenvalues[A]

y0 = 518;
r0 = 15.64;

riesenie = 
      DSolve[{y'[t] == k1*ISnerovnovaha, r'[t] == 
            k2*LMnerovnovaha, y[0] == y0, r[0] == r0}, {y, r}, 
        t]

trajektoria = 
      ParametricPlot[{y[t], r[t]} /. riesenie, {t, 0, 5}, 
        PlotRange -> {{500, 625}, {10, 16}}]

Print["t    y[t]    r[t]"]
Do[Print[
        t, ":  ", (y[t]) /. riesenie[[1]] // N, {y, r} // 
          N, "   ", (r[t]) /. riesenie[[1]] // N, {y, r} // 
          N], {t, 0, 5}]

\(<< graphics\Plotfield.m\)

g1 = 
      PlotVectorField[{Simplify[k1*ISnerovnovaha], 
          Simplify[k2*LMnerovnovaha]}, {y[t], 500, 625}, {r[t], 10, 16}, 
        PlotPoints -> 11 , ScaleFactor -> 3, 
        AspectRatio -> 0.8]

Show[g1, trajektoria]