How to write the equation of a field line of an electrostatic field? How can we write the equations of a line of force between two charges, say $q$ and $q'$?
As an example, you may consider the simpler case of two opposite charges $+q$ and $-q$, and focus on the field line emerging out of $+q$ by making an angle of $\theta$ and reaching $-q$ with an angle $\phi$ with respect to the $x$-axis (see picture below).

 A: We can find the equation for the line which makes theta angle with positive charge through these steps. I have tried to make each step clear.

A: The electric field lines are defined as being tangent in every point to the electric field in that point.
Therefore, calling $\boldsymbol r(s)$ the "trajectory" of a field line, with $s$ a parameter telling us at which point of the line we are, $\boldsymbol r(s)$ simply follows the equation
$$ \frac{d \boldsymbol r(s)}{d s} = \boldsymbol E(\boldsymbol r(s)). \tag 1$$
In your example case the electric field is given by
$$
\boldsymbol E(\boldsymbol r) = \frac{q}{4\pi \epsilon_0} \left[
\left(
\frac{x}{\left[ x^2 + y^2 \right]^{3/2}}
+ \frac{R-x}{\left[ (x-R)^2 + y^2 \right]^{3/2}}
\right) \hat{\boldsymbol x} \\
+ \left(
\frac{y}{\left[ x^2 + y^2 \right]^{3/2}}
- \frac{y}{\left[ (x-R)^2 + y^2 \right]^{3/2}}
\right) \hat{\boldsymbol y}
\right],
$$
making the solution of the system of differential equations (1) quite non trivial even in this simple case. I for one am not sure if this can be solved analytically (I had Mathematica have a try but to no avail).
If are interested in numerically verifying that this equation is true and see how the actual curve looks like, and know how to use Wolfram Mathematica, you can try the following code:
Manipulate[
 With[{
   sol = NDSolve[
     {
      x'[s] ==
       A (x[s]/(x[s]^2 + y[s]^2)^(3/2) + (
            R - x[s])/((x[s] - R)^2 + y[s]^2)^(3/2)),
      y'[s] ==
       A (y[s]/(x[s]^2 + y[s]^2)^(3/2) -
          y[s]/((x[s] - R)^2 + y[s]^2)^(3/2)),
      x[0] == 0.01,
      y[0] == 0.01 Tan[\[Theta]],
      WhenEvent[
       Abs[x'[s]] > 10^6, "StopIntegration"
       ]
      },
     {x, y}, {s, 0, 20}
     ]
   },
  ParametricPlot[
   {x[s], y[s]} /. sol,
   {s, 0, sol[[1, 1, 2, 1, 1, 2]]},
   PlotRange -> {{0, 2}, {-1, 1}}
   ]
  ],
 {{A, 0.1}, 0.001, 1, 0.01, Appearance -> "Labeled"},
 {{R, 2}, 0.001, 4, 0.001, Appearance -> "Labeled"},
 {{\[Theta], Pi/4}, -Pi/2, Pi/2, 0.001, Appearance -> "Labeled"}
 ]


A: Let a charge $+q$ be at the point $(a, 0)$ and a charge $-q$ be at the point $(-a, 0)$. Then the electric field at a point
$(x, y)$ is
\begin{equation}\tag{e1}\label{e1}
\vec{E} = q\vec{r}\left(\frac{1}{r_1^3} - \frac{1}{r_2^3}\right) - qa\hat{e}_x\left(\frac{1}{r_1^3} + \frac{1}{r_2^3}\right),
\end{equation}
where $\vec{r} = x\hat{e}_x + y\hat{e}_y$, $\vec{r}_1 = \vec{r} + a\hat{e}_x$ and $\vec{r}_2 = \vec{r} - a\hat{e}_x$. The
equation of lines of force is 
\begin{equation}\tag{e2}\label{e2}
\frac{\partial y}{\partial x} = \frac{y}{x + a\frac{r_2^3 + r_1^3}{r_2^3 - r_1^3}}.
\end{equation}
We will now solve equation \eqref{e2}. We first rearrange it as
\begin{equation}
\left((r_2^3 - r_1^3)x + a(r_2^3 + r_1^3)\right)\frac{\partial y}{\partial x} = (r_2^3 - r_1^3)y.
\end{equation}
Multiplying both sides by $y/(r_1^3 r_2^3)$, 
\begin{equation}
\left(\frac{x + a}{r_1^3} - \frac{x - a}{r_2^3}\right)y\frac{\partial y}{\partial x} = \frac{y^2}{r_1^3} - \frac{y^2}{r_2^3}.
\end{equation}
We substitute $y^2 = r_1^2 - (x + a)^2$ in the first factor on the right hand side and $y^2 = r_2^2 - (x - a)^2$ in
the second factor to get
\begin{equation}
\left(\frac{x + a}{r_1^3} - \frac{x - a}{r_2^3}\right)y\frac{\partial y}{\partial x} = \frac{1}{r_1} - \frac{(x+a)^2}{r_1^3} - \frac{1}{r_2} + \frac{(x-a)^2}{r_2^3}
\end{equation}
or,
\begin{equation}\tag{e3}\label{e3}
\frac{1}{r_1} - \frac{1}{r_2} - \frac{(x+a)}{r_1^3}\left((x + a) + y\frac{\partial y}{\partial x}\right) + \frac{(x-a)}{r_2^3}\left((x - a) + y\frac{\partial y}{\partial x}\right) = 0.
\end{equation}
We use the derivatives of $r_1$ and $r_2$ with respect to $x$
\begin{eqnarray*}
\frac{\partial r_1}{\partial x} &=& \frac{x + a}{r_1} + \frac{y}{r_1}\frac{\partial y}{\partial x} \\
\frac{\partial r_2}{\partial x} &=& \frac{x - a}{r_2} + \frac{y}{r_2}\frac{\partial y}{\partial x}
\end{eqnarray*}
in equation \eqref{e3} to get
\begin{equation}
\frac{1}{r_1} - \frac{x+a}{r_1^2}\frac{\partial r_1}{\partial x} - \frac{1}{r_2} + \frac{x-a}{r_2}\frac{\partial r_2}{\partial x} = 0,
\end{equation}
or
\begin{equation}
\frac{d}{dx}\left(\frac{x + a}{r_1}\right) - \frac{d}{dx}\left(\frac{x + a}{r_1}\right) = 0,
\end{equation}
from which we readily get
\begin{equation}\tag{e4}\label{e4}
\frac{x + a}{r_1} - \frac{x - a}{r_2} = C,
\end{equation}
where $C$ is a constant, as the solution of the differential equation (e2). This is also the solution given in article 63 of 'The Mathematical Theory of Electricity and Magnetism' by Sir James Jeans (5th edition).
A: I'll give you the start.  You need to finish the work.  
Start with an understanding of how the lines are constructed.  The lines are paths through a vector field that are resultant electric force from the particles (number does not matter) in the system.  You then connect between the charges by starting a small distance from one of the particles and move a small amount in the direction of the electric field and then take successive steps.  Eventually you will arrive at either your calculation boundary or another particle.  The equations can be derived in this mathematical method but there are going to be some errors that will vary with the size of the step you are taking.
But if you have some more math skills you can calculate the vector field and calculate the differential geometric space associated with it.  Then the equations can be directly calculated using the geodesic equation through that curved space.  That method would yield a more accurate result but it would be fun to calculate.
