How to write the equation of a field line of an electrostatic field?

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.

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}{4pi 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"}
 ]

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.

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} = qvec{r}left(frac{1}{r_1^3} - frac{1}{r_2^3}right) - qahat{e}_xleft(frac{1}{r_1^3} + frac{1}{r_2^3}right), end{equation} where $vec{r} = xhat{e}_x + yhat{e}_y$, $vec{r}_1 = vec{r} + ahat{e}_x$ and $vec{r}_2 = vec{r} - ahat{e}_x$. The equation of lines of force is begin{equation}tag{e2}label{e2} frac{partial y}{partial x} = frac{y}{x + afrac{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)yfrac{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)yfrac{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) + yfrac{partial y}{partial x}right) + frac{(x-a)}{r_2^3}left((x - a) + yfrac{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).