a = {M - kf P/(.1+P+2P^2) - kD P - kL HeavisideTheta[sin[theta]] P, eps (1/(1+P^4) - M), Pi/12}
b = {P, M, theta}
D[a,{b}]

TeXForm[%]

a = {x1^3 + 2 x2^2, 3 x1^4 + 7 x2}
b = {x1, x2}
D[a,{b}]

f = {Sin[x y] - y,y+x}

Sin[1 2]

f{2,2}

StreamPlot[f,{x,-Pi,Pi},{y,-Pi,Pi},
StreamPoints -> {{{{-2, 1}, Green}, {{0.5, -1}, Red},Automatic}}]

NTauto = {2.1(y - x/(0.1+x+2 x^2) - 0.05 x),2.1*0.05(1/(1+x^4) - y)}

StreamPlot[NTauto,{x,0,5},{y,0,1},VectorPoints->Automatic]

a = {2.1(M - P/(.1+P+2P^2) - 0.05 P), 0.105 (1/(1+P^4) - M)}
b = {P, M}
F = D[a,{b}]/. P->1.167556

TeXForm[%]

Eigenvalues[F]

f[x_, y_] = 2.1(y - x/(0.1+x+2 x^2) - 0.05 x)
g[x_, y_] =2.1*0.05(1/(1+x^4) - y)

D[f[x,y], {{x, y}}] /. x -> 1.167556

D[g[x,y], {{x, y}}] /. x -> 1.167556

vp = VectorPlot[{f[x, y], g[x, y]}, {x, 0, 4}, {y, 0, 1},
   VectorScale -> {0.045, 0.9, None},
   Axes -> True, AxesLabel -> {x, y},
   VectorPoints -> 16];
cp = ContourPlot[{f[x, y] == 0, g[x, y] == 0}, {x, 0, 4}, {y, 0, 1}];
ptRules = NSolve[{f[x, y] == 0, g[x, y] == 0}, {x, y}];
Show[vp, cp, 
 Graphics[{Red, PointSize[Large], Point[{x, y}] /. ptRules}]]

vp = StreamPlot[{f[x, y], g[x, y]}, {x, 0, 4}, {y, 0, 1},
Axes->True,AxesLabel->{x,y},StreamPoints->32];
cp = ContourPlot[{f[x, y] == 0, g[x, y] == 0}, {x, 0, 4}, {y, 0, 1}];
ptRules = NSolve[{f[x, y] == 0, g[x, y] == 0}, {x, y}];
Show[vp, cp, 
 Graphics[{Red, PointSize[Large], Point[{x, y}] /. ptRules}]]

Show[%36,AxesStyle->Black]

Show[%37,AxesStyle->Gray]

f = 2.1(y[t] - x[t]/(0.1+x[t]+2 x[t]^2) - 0.05 x[t])
g =2.1*0.05(1/(1+x[t]^4) - y[t])
F = {f,g}
DF= D[F,{{x[t],y[t]}}]

eqn[init_List] = {x'[t] == f, y'[t] == g,x[0] ==init[[1]], y[0] == init[[2]]};

eta = {1,1};
soln = NDSolve[eqn[eta],{x[t],y[t]},{t,0,240}];
ParametricPlot[Evaluate[{x[t], y[t]} /. soln], {t, 0, 240}]
Plot[Evaluate[{x[t],y[t]} /. soln], {t, 0, 240},
PlotStyle -> Automatic]