Contour Plot of system of differential equation

I am trying to find the fixed points and plot the nullcline/contourplot for the system below, but I get conditional cases

$$\frac{d\theta_1}{dt}= K sin(\theta_1 - \theta_2) - sin(\theta_1)$$ $$\frac{d\theta_2}{dt}= K sin(\theta_2 - \theta_1) - sin(\theta_2)$$

I get conditional cases when I try to use NSolve[] to find the fixed points and a very weird contour plot that I don't expect. I know that the system goes under a super critical pitchfork bifurcation at $$K = \frac {1}{2}$$ I noticed this after I was able to plot the phase portrait of the system

So, there are 3 fixed points in total.

Any help in finding the fixed points and plotting the contour plot will be appreciated. Below is my attempt

Solve[{Sin[x - y] - Sin[x] == 0, Sin[y - x] - Sin[y] == 0}, {x, y}]
ContourPlot[{Sin[x - y] - Sin[x] == 0,
Sin[y - x] - Sin[y] == 0}, {x, -4, 4}, {y, -4, 4}]
f[x_, y_, k_] := k Sin[x - y] - Sin[x];
g[x_, y_, k_] := k Sin[y - x] - Sin[y];
StreamPlot[{f[x, y, .1],
g[x, y, .1]}, {x, -\[Pi], \[Pi]}, {y, -\[Pi], \[Pi]},
VectorScale -> {0.045, 0.9, None}, Axes -> True, AxesLabel -> {x,
y}]
StreamPlot[{f[x, y, 1],
g[x, y, 1]}, {x, -\[Pi], \[Pi]}, {y, -\[Pi], \[Pi]},
VectorScale -> {0.045, 0.9, None}, Axes -> True, AxesLabel -> {x,
y}]


• There are infinitely many equilibria, which Solve returns in parametrized form, with parameter C[1] allowed to range over the integers. To get a finite number, put constraints on x and y to define a finite domain (e.g., 0 <= x < 2 Pi and so forth). Commented Jul 27, 2019 at 21:50
• The "nullclines" of the ContourPlot correspond to the vertical and horizontal tangents to the stream lines in the last StreamPlot. They are consistent, even if they are not what your expected. Commented Jul 27, 2019 at 21:53

As I mentioned in a comment, there are infinitely many equilibria, which Solve[] returns in parametrized form, with parameter C[1] allowed to range over the integers. To get a finite number, put constraints on x and y to define a finite domain (e.g., -Pi < x < Pi and so forth to accord with the plots). We can put the elements together in a single plot that makes a nice illustration of the flow of the phase field.

cs = {Directive[Dashed, AbsoluteThickness[1], Hue[0.9]],  (* contour styles *)
Directive[Dashed, AbsoluteThickness[1], Hue[0.1]]};
legend = Style[LineLegend[Append[cs, Red],                (* legend *)
{RawBoxes@RowBox[{"\[Piecewise]",
GridBox[{{ToBoxes[\[DifferentialD]x == 0]}, {"Vertical Tangent"}}]}],
RawBoxes@RowBox[{"\[Piecewise]",
GridBox[{{ToBoxes[\[DifferentialD]y == 0]}, {"Horizontal Tangent"}}]}],
RawBoxes@RowBox[{"\[Piecewise]",
GridBox[{{ToBoxes[\[DifferentialD]x == \[DifferentialD]y == 0]}, {"Equilibria"}}]}]
},
Joined -> {True, True, False}], GridBoxOptions -> {ColumnAlignments -> Left}];

f[x_, y_, k_] := k Sin[x - y] - Sin[x];
g[x_, y_, k_] := k Sin[y - x] - Sin[y];
pts = Solve[{f[x, y, 1] == 0, g[x, y, 1] == 0, -π < x < π, -π < y < π}, {x, y}];
nc = ContourPlot[{f[x, y, 1] == 0, g[x, y, 1] == 0},
{x, -π, π}, {y, -π, π}, ContourStyle -> cs];
sp = StreamPlot[{f[x, y, 1], g[x, y, 1]},
{x, -π, π}, {y, -π, π}, VectorScale -> {0.045, 0.9, None}];

Legended[
Show[
nc, sp, Graphics[{Red, AbsolutePointSize[6], Point[{x, y} /. pts]}],
FrameLabel -> {x, y}],
legend
]


• thanks a lo for the clarification! Commented Jul 28, 2019 at 17:12
• @Rumman You're welcome. Commented Jul 28, 2019 at 17:57