Let's use the convention for the wavefunction of the hydrogen atom as in Wikipedia. We're interested in the second-order correction to the energy of the ground state due to an electric field. In appropriate units the ground state Hamiltonian is $H_0 = -\frac{1}{2}\nabla^2 -\frac{1}{r}$. If we add an electric field in the $z$ direction, the Hamiltonian becomes $H=H_0+H_1$ with $H_1 = Fr\cos{\theta}$. One can choose $F$ positive or negative depending on whether the field is in the positive or negative $z$ direction. The sign choice we make differs from OP's unnumbered equation so that in the end we get $A$ as in the OP's question; if we chose $H_1$ as in the OP's question, then the sign in $A$ would be flipped. I will follow Aruldhas section 9.5.
We're interested in the shift in the ground state, that is the $n=1, l=0$ (and therefore $m=0$) state of the hydrogen atom. We write $\phi^{(0)}_n$ for the $n$-th $l=0$ eigenstate of $H_0$. By direct computation, the first order energy shift vanishes. The second order energy shift is given by
$$E^{(2)}_{1}=\sum_{m, m\neq 1}\frac{|\langle\phi_m^{(0)}|H_1|\phi_1^{(0)}\rangle|^2}{E_1^{(0)}-E_m^{(0)}}$$
This should be compared with OP's eq. $(1)$. However this can be written in a different form which can be more suited to our uses. The first order correction to the ground state wavefunction is
$$\phi^{(1)}_1= \sum_{m, m\neq 1}\frac{\langle\phi_m^{(0)}|H_1|\phi_1^{(0)}\rangle}{E_1^{(0)}-E_m^{(0)}}\phi_m^{(0)}$$
so that the corrected wavefunction (to first order) would be $\phi^{(0)}_1+\phi^{(1)}_1$. Therefore one can also write the second order energy shift as
$$E^{(2)}_{1}=\langle\phi_1^{(0)}|H_1|\phi_1^{(1)}\rangle$$
This formula is not very useful as it is. But now suppose that we could find an operator $A$ such that
$$A|\phi_1^{(0)}\rangle = |\phi_1^{(1)}\rangle$$
Then we get $E^{(2)}_{1}=\langle\phi_1^{(0)}|H_1 A|\phi_1^{(0)}\rangle$, which is substantially simpler than evaluating the full infinite sum, at least if $A$ is simple enough. To determine $A$ we can't simply use its definition, because we would need to compute $\phi_1^{(1)}$ and that's as bad as computing the infinite sum. So we take a different route. The eigenvalue equation $H \phi = E \phi$ can be written to first order (not zeroth order) as
$$H_1\phi_1^{(0)} +H_0\phi_1^{(1)} =E_1^{(1)} \phi_1^{(0)} +E_1^{(0)}\phi_1^{(1)}$$
But the first order energy shift $E_1^{(1)}$ vanishes, so using $A$ and $H_0$ we get
$$(AH_0-H_0A)\phi_1^{(0)} = H_1\phi_1^{(0)}\tag{2'}$$
This should be compared with OP's eq. $(2)$, the main difference being that the relation must only be satisfied for the $n=1$ state (and in fact, as stated, OP's eq. $(2)$ is not true for $n\neq 1$). Note that if $A$ is a solution of $(2')$ and $B$ is any operator which commutes with the hamiltonian, then $A+B$ will also be a solution of $(2')$. Now as a simplifying assumption, suppose that $A$ is a scalar function (i.e. it contains no derivatives). After all, we don't care about a specific solution, any solution will do. Since $A$ is a scalar, it commutes with the $1/r$ part of $H_0$. So expanding $(2')$ using the chain rule we get
$$\phi_1^{(0)}\nabla^2 A +2 \nabla A \cdot \nabla \phi_1^{(0)}=2H_1 \phi_1^{(0)}$$
In OP's choice of units and with the Wikipedia convention, $\phi_1^{(0)} = C e^{-r/n}$ where $C$ is some constant. So in fact we get
$$\nabla^2 A - 2\frac{\partial A}{\partial r}=2H_1$$
As for how to solve this, it is usually very hard to solve PDEs and in many cases most of it is guesswork. The only $\theta$ dependence in $H_1$ is through $\cos{\theta}$ which is equal to $P_1(\cos{\theta})$, where $P_1$ is a Legendre polynomial. These are eigenfunctions of the laplacian, so we can guess $A = B(r) \cos{\theta}$ and in fact the $\theta$ dependence vanishes, and we get
$$\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial B}{\partial r}\right)-2\frac{1}{r^2}B-2\frac{\partial B}{\partial r}=2Fr$$
or expanding
$$\frac{\partial^2 B}{\partial r^2}+2\left(\frac{1}{r}-1\right)\frac{\partial B}{\partial r}-2\frac{1}{r^2}B=2Fr$$
which is an equation for $B$ only. We can then look for a series solution $B(r)=F \sum_l a_l r^l$, which gives $A$ as in the OP's question. Alternatively, we see that the coefficients of the equation are singular at $r=0$. Therefore expanding for small $r$ we get the equation
$$\frac{\partial^2 B}{\partial r^2}+2\frac{1}{r}\frac{\partial B}{\partial r}-2\frac{1}{r^2}B=0$$
This is an equidimensional equation which has a solution of the form $r^k$. Substituting and solving the quadratic we get $k=1$ or $k=-2$. We choose $k=1$, which is non singular, so in fact $B\approx r$ near zero. It is therefore useful to write $B(r) = r F G(r)$ (we put in $F$ to get rid of it). Substituting we get a simpler equation
$$\frac{\partial^2 G}{\partial r^2}+(4-2r)\frac{\partial G}{\partial r}-2G=2r$$
This you can plug into WolframAlpha.They give the solution in a complicated form, but some things cancel out and you can choose the constants to get rid of the nasty part, getting the particular solution $G = -(r/2+1)$.
