Inv = i0 - i1*r[t];
S = (1 - b)*(y[t] - T) - a;
IS = (Inv + G == S + T);
ISN = Inv + G - (S + T);
L = c0 + c1*y[t] - c2*r[t];
LM = (L == M);
LMN = (L - M);
y'[t] == 
      Collect[Expand[k1*ISN], {y[t], r[t]}]
r'[t] == Collect[Expand[k2*LMN], {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 = 1;
k2 = 0.5;
y'[t] == Collect[Expand[k1*ISN], {y[t], r[t]}]
r'[t] == Collect[Expand[k2*LMN], {y[t], r[t]}]

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

dame pociatocnu hodnotu 518 a 15.64. a zistime, ci to konverguje ku 599.761 a 13.3002

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

    MatrixForm[A]

Eigenvalues[A]

y0 = 518;
r0 = 15.64;
riesenie = 
      DSolve[{y'[t] == k1*ISN, r'[t] == k2*LMN, 
          y[0] == y0, r[0] == r0}, {y, r}, t]

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

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

\(<< graphics\Plotfield.m\)

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

Show[g1, trajektoria]