Deriving solution of the Renormalization Group Equation

I am trying to follow Matthew Schwartz's renormalization group lectures (pdf or see Chapter 23 of QFT and the SM by Matthew Schwartz), but I am having trouble with Eq. (book 23.31/pdf 29).

I understand that the $\beta$-function comes from the renormalization group equation (RGE), $$\beta(\alpha)\equiv \mu\frac{\partial\alpha}{\partial\mu} \tag{book 23.30/pdf 28},$$ but then the notes introduce $$\beta(\alpha) = -2\alpha\biggl[\frac{\epsilon}{2} + \biggl(\frac{\alpha}{4\pi}\biggr)\beta_0 + \biggl(\frac{\alpha}{4\pi}\biggr)^2\beta_1 + \biggl(\frac{\alpha}{4\pi}\biggr)^3\beta_2 + \cdots\biggr] \tag{23.31/29},$$ where I've corrected a typo in the pdf. I totally did not understand this expansion. Then it finally goes to $$\alpha(\mu)=\frac{2\pi}{\beta_0}\frac{1}{\ln(\frac{\mu}{\Lambda_\text{QED}})}\tag{23.32/30}$$ where $\Lambda_\text{QED}$ is the Landau pole and the fundamental scale of QED. Where does it come from? I'd appreciate it if anyone could help me to derive Eq. 23.31/29 and 23.32/30.

He is writing the $\beta$ function as a perturbation series in the small parameter $\alpha$, $$\beta(\alpha) = \sum_{n=1}^\infty c_n \alpha^n.$$ By convention, it's common to write this expansion in the form he has given. The first coefficient is always $c_1 = - 4\pi\epsilon$ (where $\epsilon$ is the dimensional regularization parameter). The remaining terms define the coefficients $\beta_n$, which one computes in perturbation theory. For example, he has computed the series up to one loop for QED, which he writes in eqn. 27: $$\mu\frac{\mathrm{d}e}{\mathrm{d}\mu} = -\frac{\epsilon}{2} e+ \frac{e^3}{12\pi^2}+\cdots$$ Or, in terms of $\alpha :=e^2/4\pi$, $$\mu \frac{\mathrm{d}\alpha}{\mathrm{d}\mu} = \frac{e}{2\pi} \mu \frac{\mathrm{d}e}{\mathrm{d}\mu} = -\frac{\epsilon}{4\pi}e^2+\frac{e^4}{24\pi^3} +\cdots = -\epsilon\alpha+\frac{2}{3\pi}\alpha^2+\cdots,$$ i.e. $\beta_0 = -4/3$. He sends the regulator $\epsilon$ to 0 and solves the differential equation $$\int_{\alpha(\mu_0)}^{\alpha(\mu)}\frac{\mathrm{d}\alpha}{\alpha^2} = \frac{2}{3\pi} \int_{\mu_0}^\mu\frac{\mathrm{d}\mu}{\mu},$$ which gives $$\alpha(\mu) = -\frac{3\pi}{2}\frac{1}{\log(\mu/\mu_0)}.$$ $\mu_0$ is an integration constant which he has chosen to be $\Lambda_{\text{QED}}$.
Schwartz is simply noting that the $\beta$-function has a generic expansion in QED of the form (29) where $\beta_{0,1,2}$ are some numbers that can be computed by explicitly calculating the various Feynman diagrams.
For instance the leading $\epsilon/2$ is the tree-level result in $d=4-\epsilon$ dimensions. This can be easily seen as follows. In $d=4-\epsilon$, we have $[e]=\epsilon/2$ (do you know how to see this??). Then, we find $[\alpha] = [e^2] = \epsilon$. Thus, we can define the dimensionless renormalized coupling $$\alpha(\mu) = \mu^{-\epsilon} \alpha$$ This is of course a tree-level result as generically $\alpha$ obtains an anomalous dimension. From this equation, we can determine $$\beta(\alpha) = \mu \frac{d\alpha}{d\mu} = - \epsilon \alpha$$ This is precisely the first term that he has in the $\beta$-function expansion. If I then take into consideration quantum renormalization effects the $\beta$-function takes the more general form (29).
Good! Now, we understand (29). From (29), we can derive (30). Keeping terms to leading order in $\alpha$ in $d=4\implies \epsilon=0$, we find $$\beta (\alpha) = \mu \frac{d\alpha}{d\mu} = - \frac{\alpha^2}{2\pi} \beta_0$$ This differential equation can be solved. Doing so introduces an integration constant $\Lambda_{QED}$ and the solution is of the form $$\alpha(\mu) = \frac{2\pi}{\beta_0} \frac{1}{\log \frac{\mu}{\Lambda_{QED}}}$$ We can now interpret $\Lambda_{QED}$ as the energy scale where $\alpha(\mu) \to \infty$, i.e. it is the scale where the UV description of the theory breaks down. It is known as the Landau pole.