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 = 0;
c1 = 0.25;
c2 = 1;
M = 450;
i0 = 0;
i1 = 3;
a = 400;
b = 0.8;
T = 0;
G = 0;
k1 = 0.4;
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 = 1900;
r0 = 10;
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, 100}, 
        PlotRange -> {{1830, 1900}, {8, 19}}]

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], 1830, 
          1900}, {r[t], 8, 19}, PlotPoints -> 11, ScaleFactor -> 3, 
        AspectRatio -> 0.8]

Show[g1, trajektoria]