solution of $y' + y^2 = \varphi^2(x)$ I need to solve differential equation in the interval $[-\pi/2,\pi/2]$
\begin{eqnarray}
y''(x) = y(x)\sin^2x
\end{eqnarray}
Trying $y(x) = \exp(\psi(x))$  yields,
\begin{eqnarray}
\zeta'(x) + \zeta^2(x) = \sin^2x \hspace{2cm} \zeta(x) = \psi'(x)
\end{eqnarray}
This equation seems to simpler than the original second order equation but
still I can't find way to solve this. Even if the equation is solved for function other than $\sin(x)$ with some important properties of $\sin(x)$ preserved,  I will consider myself fortunate.
Is it possible to solve for any $\varphi(x)$ such that,
\begin{eqnarray}
\zeta'(x) + \zeta^2(x)  = \varphi^2(x)
\end{eqnarray}
Where, $\varphi(x)$ is a monotonic function some interval $[a,b]$, with exactly one inflation point at $(a+b)/2$ and derivative vanishing at endpoints?
\begin{eqnarray}
\varphi'(x) >= 0 \\
\varphi''(x)|_{x=\frac{a+b}{2}} = 0\\
\varphi'(x)|_{x=a,b} = 0
\end{eqnarray}
 A: As Andrew D pointed out, the Mathieu differential equation is:
$$y''(x)+\left(a-2q\cos (2x)\right)y=0$$
When $a=-1/2$ and $q=-1/4$, this gives your original equation.
Thus, given the complex valued Mathieu functions $C(a,q,x), S(a,q,x)$:
$$y(x)=k_1 C(-1/2,-1/4,x)+k_2 S(-1/2,-1/4,x)$$
Where $k_1 = y(0) / C(-1/2, -1/4,0)$ and $k_2 = y'(0) / S'(-1/2, -1/4,0)$, and $S'$ is the derivative of the Mathieu $S$ function. You can show that $y(x)$ is real for all real $y(0),y'(0)$.
A: Let 
$$
v = y’
$$
then
$$
y’’ = \frac{dv}{dx} = \frac{dy}{dx} \frac{dv}{dy} = v \frac{dv}{dy}
\\ \frac{y''}{y} = \frac{v}{y} \frac{dv}{dy} = \frac{dv^2}{dy^2} = \sin^2{x}
$$
This might be used as the basis of a series solution. For example let
$$
v = (a_0 + a_1y + a_2y^2 + ...)^{-1}
$$
Then the LHS becomes
$$
\frac{dv^2}{dy^2} = \frac{v}{y} \frac{dv}{dy} =-\frac{(a_1 + 2a_2y + 3a_3y^2...)}{y(a_0 + a_1y + a_2y^2 + ...)^3}
$$
and
$$
v = \frac{dy}{dx}= (a_0 + a_1y + a_2y^2 + ...)^{-1}
$$
can be integrated to obtain
$$
x = a_0y + \tfrac{1}{2}a_1y^2 + \tfrac{1}{3}a_2y^3 + ... 
$$
where the constant of integration has been set to $0$ to fix the $x$ origin where $y=0$, and $a_0$ is determined by the gradient of $y(x)$ at $x=0$:
$$
a_0 = 1/y'(0)
$$
We now have
$$
-\frac{(a_1 + 2a_2y + 3a_3y^2...)}{y(a_0 + a_1y + a_2y^2 + ...)^3} = \sin^2{(
a_0y + \tfrac{1}{2}a_1y^2 + \tfrac{1}{3}a_2y^3 + ...)}
$$
Expressing the RHS as a Taylor series and equating coefficients of $y^n$:
$$
a_1 =0, a_2 = 0, a_3= 0, a_4=-\tfrac{1}{4} a_0^5
$$
and so on.
So as a first approximation
$$
x \approx a_0y + \tfrac{1}{5}a_4y^5 = a_0y - \tfrac{1}{20}a_0^5y^5 
$$
