Consider the Schrödinger equation: $$y''(x)+EQ(x)y(x)=0,\tag{1}$$ where $E>0, Q(x)>0, y(0)=y(\pi)=0. $
Use WKB approximation to obtain $$y(x) \sim CQ^\frac{-1}{4}(x)\sin\Big(\sqrt(E) \int_{0}^{x} \sqrt(Q(t) dt \Big),\ E \to \infty. $$
What I did:
Let $\phi(x;\epsilon)=\log (y(x;\epsilon)).$ So, $y(x;\epsilon)=exp(\phi(x;\epsilon).$
Substitute this to the main equation to get: $$\phi''(x;\epsilon)+\phi'(x;\epsilon)^2+EQ(x)=0. \tag{2}$$
Suppose $$\phi(x;\epsilon)=\frac{1}{\epsilon}\sum_{n=0}^{\infty} \phi_n(x)\epsilon^{n},\ \epsilon \to 0. $$
Substitute this to $(2)$ using $\phi(x;\epsilon)=\phi_0(x) \frac{1}{\epsilon}+O(1)$ we get,
$$\frac{1}{\epsilon}\phi''_0(x)+\frac{1}{\epsilon^2}\phi'_0(x)^2+EQ(x)+O(1)=0.$$
This implies $$\epsilon \phi''_0(x)+\phi'_0(x)^2+\epsilon^2EQ(x)+O(\epsilon^2)=0.$$
So, at $O(1)$, $$\phi'_0(x)^2=0.$$
This means $\phi_0(x)=C$ where $C$ is constant. We then continue to the next order by substituting $\phi(x;\epsilon)=\phi_0(x) \frac{1}{\epsilon}+\phi_1(x)+O(\epsilon)=\frac{C}{\epsilon}+\phi_1(x)+O(\epsilon)$ to $(2)$.
We get $$\phi''_1(x)+\phi'_1(x)^2+EQ(x)=0.$$
This is where I'm stuck, any idea how to solve this ODE to get $\phi_1(x)$?
No, your first balancing decision does not make sense. In the balancing you attempt to get two terms of equal magnitude/scale. You only got one dominating term, which you saw gave a trivial result in $ϕ_0=const.$ and no simplification relative to the original equation in the next equation for $ϕ_1$.
You need to include in the balancing calculations that $E$ is assumed large, and $ϵ$ small so that the balancing of the largest terms among the scale coefficients $ϵ,1,ϵ^2E$ in $$ϵ[ϵϕ''(x;ϵ)]+[ϵϕ'(x;ϵ)]^2+ϵ^2E[Q(x)]=0$$ gives $ϵ^2E=1$. Then the equation for the first term $ϕ_0$ in the expansion of $ϕ$ is non-trivial $$ ϕ_0'(x)^2+Q(x)=0. $$ The next order equation gives $$ ϕ_0''(x)+2ϕ_0'(x)ϕ_1'(x)=0\implies ϕ_1'(x)=-\frac{ϕ_0''(x)}{2ϕ_0'(x)}= -\frac{Q'(x)}{4Q(x)}=-\frac14(\ln|Q(x)|)' $$ which when integrated gives the cited WKB approximation formula.
The next order term has the equation $$ ϕ_1''(x)+ϕ_1'(x)^2+2ϕ_0'(x)ϕ_2'(x)=0\implies ϕ_2'(x) = \frac{ϕ_0'''(x)}{4ϕ_0'(x)^2} - \frac{3ϕ_0''(x)^2}{8ϕ_0'(x)^3}=\left(\frac{ϕ_0''(x)}{4ϕ_0'(x)^2}\right)'+\frac{ϕ_0''(x)^2}{8ϕ_0'(x)^3} $$ which does not look like a nice-to-integrate term.