Trajectories that connect equilibrium points

Your question is very general, but I think it's possible to give recommendations that work in many cases. Although they are quite general too :)

There are two methods for solving such problems that immediately come to my mind. The first one is finding trapping region and the second one uses the existence of first integral. I should note that the first method doesn't compute the trajectory: it just shows that there's a trajectory that connects two equilibria. However, the second method might give you an implicit equation for finding separatrix.

Let's move to the details.

Finding trapping region

I hope that this simple illustration will make things clear:

The idea basically is the following. You have a compact domain $\overline{D}$ which is forward-time flow invariant and one of separatrices stays in this domain. Since domain is forward-time flow invariant, this separatrix can't escape it. According to Poincare-Bendixson theory the only limit sets that are possible for trajectories that stay in bounded domains are equilibria, limit cycles or separatrix contours with saddles. If you somehow know that there is only one pair of equilibria in $\overline{D}$ and no limit cycles then the only possible conclusion is that separatrix goes to this other equilibrium.

How to determine if domain is trapping region or not? If you have an analytical expression for the boundary of domain, you can always check the dot product of vector field and normal vector along the boundary. Depending on the sign of this product and on orientation of boundary you can tell whether vector field points inward or outward. If vector field points strictly inward at the boundary $\partial D$ of domain — then it's a trapping region. If it's pointing inward except at few equilibria (like on the picture) — it's also a trapping region. A particular example of finding such trapping region with similar picture can be found here.

Knowing first integral

Sometimes it's possible to find locally non-constant smooth function $F(x,y)$ such that $\frac{d}{dt} F(x,y) \equiv 0$. This function is called "first integral". The implication of existence of first integral is crucial for dynamics: trajectories of system lie on level sets of $F(x,y)$. Another important fact is that having the first integral makes impossible to have sources or sinks: by continuity whole sink's domain of attraction would have the same value of $F(x,y)$ as the sink itself, which would contradict assumption that $F(x,y)$ is locally non-constant. Level sets of $F(x,y)$ consist of whole trajectories and maybe of union of them with equilibria. It's easy to prove that separatrices of saddle lie in the same level set as the saddle itself (straightforward conclusion from continuity). Also, if separatrix goes to some equilibrium then it also lies in the same level as the separatrix. So this gives a sort of an algorithm: find equilibria, check if any of them lie in the same level set by computing $F(x,y)$ at them. Then just study how this level set looks like. So, you might have separatrices connecting two different saddles (heteroclinic connection) or connecting one saddle with itself (homoclinic loop).

Someone constructed system with heteroclinic trajectory in mind

The explanation for this paragraph also can be used as a method of constructing systems with heteroclinic trajectory. Suppose you have a system

$$\dot{x} = P(x, y), \; \dot{y} = Q(x, y)$$

for which you know that set $H(x, y) = 0$ contains heteroclinic or homoclinic trajectory with equilibria that are joined by it. Then the vector field defined by system

$$\dot{x} = P(x, y) + H(x, y) \cdot \tilde{P}(x, y), \; \dot{y} = Q(x, y) + H(x, y) \cdot \tilde{Q}(x, y)$$

is the same at $H(x, y) = 0$ as the vector field of previous system, hence they share phase trajectories that lie on $H(x, y) = 0$. Thus, the second system also has heteroclinic or homoclinic trajectory. This method is just perfect when the initial system has first integral $\mathcal{F}(x, y)$ which can be used for proving the existence of heteroclinic solution. You can always define $H(x, y)$ in such way that heteroclinic trajectory lies in $H(x, y) = 0$. But maybe you have multiple heteroclinic trajectories that you want to preserve? Alright then, go for $$H(x, y) = (\mathcal{F}(x, y) - C_1)\cdot \ldots \cdot (\mathcal{F}(x, y) - C_n),$$ where $C_1, \dots, C_n$ are level sets where your heteroclinic trajectories lie. This method can be applied not only for preserving heteroclinic or homoclinic trajectories, but also for preserving closed trajectories to obtain examples with limit cycle.

There also two more "esoteric" methods :)

Knowing Lyapunov function in domain that contains saddle and sink

Suppose that you have a sink for which you have constructed a Lyapunov function in some domain. Suppose that you were so lucky that separatrix enters this domain. Then everything is done — separatrix goes to sink. Why is it "esoteric"? Finding trapping region (generally) is hard, finding Lyapunov function with big domain is much harder (IMHO, of course). However, existence of both things is guaranteed for some cases in qualitative theory.

There's no first integral, but you somehow found one-dimensional invariant manifold

Sometimes you have invariant axes. Or just invariant straight lines. Sometimes you have equilibria on these lines and some of them are saddles. Well, again, you're done. (I'm speaking about straight lines because they're easiest for check) Maybe sometimes you'll be lucky and you'll find some invariant graph $(x, f(x))$ (there's a simple condition for this graph to be invariant, but it maybe not so easy to fulfill them), where you can study the restricted dynamics and find saddles again. I think you've got an idea. See how many "maybe" and "sometimes" ? Yep, that makes this method "esoteric" too.

All questions, comments and corrections are welcome.