I have a differential equation of the form
$$ \frac{\text{d}y(x)}{\text{d}x} - c_1z_1(x)x^2 \int_x^{c_2} \frac{z_1(x)}{x^2 z_2(x)} y \,\text{d}x = -z_3(x) $$
where $y, z_1, z_2$ and $z_3$ are functions of $x$ and $c_1$ and $c_2$ are constants. I have specific values of $x$ and cannot evaluate these functions at arbitrary locations, although I could interpolate. Further, $x$ is bounded between, say, 0 and $1=c_2$.
I am wondering what would be a good method for numerically solving this equation for $y(x)$?
This is not a simple differential equation, but an integro-differential equation, which are typically much more complicated to solve. The variational iteration method is usually described as simple. See Wang S.-Q. and He J.-H., Variational iteration method for solving integro-differential equations, Phys. Lett. A., 2007. Note that to numerically solve anything, you will need an initial condition $y(x_0) = a$.
Edit : As pointed out by comments, this can be written as a second order equation, in which case usual differential algorithms can be used. The simplests one are finite differences, where $y''(x) = \Delta^{-2} \frac{y(x-\Delta) - 2 y(x) + y(x+\Delta)}{2}$ and $y' = \Delta^{-1} \frac{y(x+\Delta) - y(x-\Delta)}{2}$. Solving the equation on a grid X will generate a nonlinear equation which in matrix form should look as $A(X) X = B(X)$, which you can typically solve using Newton's method.
EDIT2: A few more precisions on getting the matrix form. Assuming you have a grid of points $\vec{x} = [x_1, x_2, ..., x_{i-1}, x_i, x_{i+1}, ..., x_n]$ and $\vec{y} = [y_1, y_2, ..., y_{i-1}, y_i, y_{i+1}, ..., y_n]$, then the operator $y''(x) =1/2 \Delta^{-2} (y(x+\Delta) -2y(x) + y(x-\Delta))$ on the grid is given by $y''(x_i) = \Delta^{-2}/2 * (y_{i-1} -2 y_i + y_{i+1})$, so that the whole operator $y''(x)$ is a sparse tridiagonal matrix $M$ consisting of $diag(M) = -\Delta^{-2}$ and offdiagonals $diag_{-1}(M) = diag_{+1}(M) = \Delta^{-2}/2$. Th operator $y'(x)$ can be represented in the same way. For the grid points where you want to enforce boundary conditions, namely $y_1$ and $y_n$, your equation is different and only reads $y_1 = a_1$ and $y_n = a_n$. You can enforce this by setting the matrix coefficients to be one on the diagonal ($A_{1,1} = 1$, $A_{n,n} = 1$) as well as the RHS $B_1 = a_1$ and $B_n = a_n$. Enforcing derivatives is done in the same way by using first order approximations $y'(x_1) = \Delta^{-1} (y_2 - y_1)$, which corresponds to $A_{1,2} = \Delta^{-1}, A_{1,1} = -\Delta^{-1}$. If your system is linear in $y$ then your solution is given by a simple linear equation solving problem, as it becomes $A \vec{y} = B$ and $\vec{y} = A^{-1} B$. If you're using MATLAB, consider using A\B to get the solution.