Autonomous system of 2 non-linear differential equations Hello I am trying to solve this system:
$$ \alpha '(t)= -\frac{e^{-3 \alpha (t)} \sin (2 \mu p_\alpha(t))}{\mu} \tag{1}\label{1}$$
$${ p_\alpha}'(t)= 3 \rho _0^2 (w-1) e^{-3 w \alpha (t)} \tag{2}\label{2} $$
Where $\alpha(t)$ and $p_\alpha(t)$ are the two functions of time that I am trying to calculate. Actually I am interested only in $\alpha(t)$. $\rho_0,w,\mu$ are constant parameters where $\rho_0>0$, $-1<w<1$ and $\mu>0$ is a small perturbation.
As you might have already guessed these are nothing else than the Hamilton's equations of a dynamical system endowed with an Hamiltonian. The Hamiltonian of the system, modulo a constant, is: 
$$\mathcal{H}= e^{-3\alpha(t)}\left(p_{-}^2+p_{+}^2-\frac{\sin(\mu p_\alpha(t))^2}{\mu^2}+\rho_0^2 e^{-3(w-1)\alpha(t)}\right)=0$$
Since we are dealing with general relativity the Hamiltonian is null (frozen time formalism).
Anyway I am able to solve the system in the limit $\mu\to0$ and the result is correct (I know for sure). 
This is the procedure I follow: I divide \eqref{2} over \eqref{1} to find $\frac{d p_\alpha}{d\alpha}$. I integrate it to find $p_\alpha(\alpha)$ that I then substitute in \eqref{1}. Now I have only an ODE in one variable that I can integrate to find $\alpha(t)$.
As I said in the unperturbed case ($\mu\to0$) this procedure gives me the correct result but for $\mu\neq0$ the very same procedure gives me nonsense.
I am using Wolfram Mathematica to do the calculations, but I suspect the problem is not the code that I tried to manipulate in every conceivable (to me at least) way. 
Every kind of advice is welcomed. I don't expect you to work out the solution for me. I only would like to know how would you tackle the problem. Thank you 
 A: This is a part of the chat discussion that is too long to be put in chat. And chat doesn't render math.
In theory I might be doing mistake in my calculations, but if I hadn't then the qualitative conclusion about this system is the same as yours about cyclic change of Universe's volume. Actually, all trajectories of this system are periodic. I would be very-very glad if someone could re-check my calculations!
In principle this system could be solved, but I've decided to study its Poincaré map. It actually requires less calculations and these calculations are very simple. Because the system is periodic in $p_{\alpha}$ the phase space is a cylinder actually. This is due to the fact that the system posesses a discrete translation symmetry: if $(\xi(t), \chi(t))$ is a solution then $(\xi(t), \chi(t) + 2\pi k)$ is also a solution for any $k \in \mathbb{Z}$. This system has no equilibria so its trajectories are likely to wind up around the cylinder. Note that if you take a line $p_{\alpha} = \frac{\pi k}{\mu}$ the vector field is strictly orthogonal to it — it's one of the nullclines and a good candidate for a Poincaré section. Thus we can represent cylinder as a strip $- \frac{\pi}{\mu} \leqslant p_{\alpha} \leqslant 0$ with glued lines $ p_{\alpha} = 0$ and $p_{\alpha} = -\frac{\pi}{\mu}$. 
The reason why I've chosen exactly this strip is the following. Note that  $p'_{\alpha} < 0$ everywhere because of constraints on values of parameters. If we start at $p_{\alpha} = 0$ we go below until we hit $p_{\alpha} = -\frac{\pi}{\mu}$. Basically this would be our Poincaré map $\mathcal{F}$: start from $(\alpha, 0)$, go with the flow until you hit $p_{\alpha} = -\frac{\pi}{\mu}$ at point $(\mathcal{F}(\alpha), -\frac{\pi}{\mu})$. The good thing is that in your case the Poincaré map can be explicitly computed.
How will we do it? We won't write the solution of the system, solve for arrival time and plug it into the solution. Instead we will use equation 
$$ \frac{d \alpha}{d p_{\alpha}}= -\frac{e^{-3 \alpha} \sin (2 \mu p_\alpha)}{3\mu \rho^2_0 (w-1) e^{-3w\alpha}}.  $$
What is it good for? Why we can use it? Well, $t$ changes strictly monotonously, $p_{\alpha}$ changes in the same fashion (and derivative is never zero). By some sort of inverse theorem we can now consider $p_{\alpha}$ as our new independent variable. This is very convenient: the equation is separable and Poincaré map can be written simply as a solution of some Cauchy problem. Let's do it: if trajectory starts at $(u, 0)$ and goes to $(\mathcal{F}(u), -\frac{\pi}{\mu})$, then the solution is 
$$ -3(w-1) e^{-3\alpha(w-1)} \,d\alpha = \frac{\sin{(2\mu p_{\alpha})}}{\mu\rho^2_0}\, d p_{\alpha}$$
$$ \int\limits_{u}^{\mathcal{F}(u)}-3(w-1) e^{-3\alpha(w-1)} \,d\alpha = \int\limits_{0}^{-\frac{\pi}{\mu}}\frac{\sin{(2\mu p_{\alpha})}}{\mu\rho^2_0}\, d p_{\alpha}$$
$$\left . \left (e^{-3\alpha (w-1)} \right ) \right \vert^{\mathcal{F}(u)}_{u} = \left . \left ( - \frac{\cos{2 \mu p_{\alpha}}}{2\rho^2_0 \mu^2} \right ) \right \vert^{-\frac{\pi}{\mu}}_{0}$$
The function in the RHS is $\frac{\pi}{\mu}$-periodic and we take the integral over the period so it's zero (the same conclusion follows from direct substitution) and we arrive at equation
$$e^{-3\mathcal{F}(u) (w-1)} - e^{-3{u} (w-1)}   = 0 $$
From this we get that $\mathcal{F}(u) = u$ which means that every solution is periodic. It doesn't necessarily mean that solution $\alpha(p_{\alpha})$ is constant — no, we can check this directly and it's not so. All solutions are just some closed curves that wrap around the cylinder. They don't necessarily have the same period w.r.t. to time $t$ (although they have the same period w.r.t. $p_{\alpha}$ when it's considered as independent variable). But alas, this map is an identity and doesn't have a good non-trivial dynamics which could be in accordance with your colleague's observations.
