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

$\begingroup$

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"}
 ]