I was trying to derive the Bernoulli equation from the above equation
for time independent flow
If you are studying something time independent then you just let $\frac{\partial}{\partial t}$ to be zero:
$$
\frac{\partial}{\partial x_j} \left [ \frac{1}{2} \rho v^2 v_j + \rho h v_j + \rho \phi v_j \right ]=0
$$
Next step is to get rid of $\rho$. Bernoulli's equation doesn't contain $\rho$, does it? We'll need mass balance for the stationary state:
$$\frac{\partial \rho v_j}{\partial x_j} = 0$$
which together with the first equation leads to:
$$
\rho v_j \frac{\partial}{\partial x_j} \left [ \frac{1}{2} v^2 + h + \phi \right ]=0
$$
Now one should recall that Bernoulli's law is valid only along streamlines/pathlines, which do coincide for a steady flow. If something called $A$ is constant along the vector field $\boldsymbol v$, then it should suffice the equation
$$v_j \frac{\partial A}{\partial x_j} = 0$$
Actually it is just a derivative of $A$ along the vector field $\boldsymbol v$ --- if the derivative is zero then $A$ is constant along the vector field.
Thus we got the Bernoulli's law:
$$
\frac{1}{2} v^2 + h + \phi = const
$$
along the streamline/pathline for the stationary case.
If the contents of the whole pipe are exposed to the atmosphere then by definition there is no pressure differential between inlet and outlet except caused by viscous drag. I am assuming you are disregarding this considering your use of the Bernoulli equation. Since $P(t)=0$ it follows quite simply that the answer to your question is the solution to the Ricatti equation:
$$\frac{dv(t)}{dt}+av^2=c$$
If $P(t)\neq0$, then the form of $P(t)$ needs to be known otherwise the equation is simply not solvable.
The question arises if the (unsteady) Bernoulli equation is at all applicable in this case, the decreasing velocity profile in the pipe is mostly determined by the viscous drag in the fluid and at the walls. Considering viscous effects, Bernoulli's equation is not applicable unless an additional term is added to account for this:
$$\frac{d}{dt}(Me)=-\Delta\phi_{m}e+\phi_w-\epsilon$$
where $e=\frac{p}{\rho}+\frac{1}{2}v^2+gz$ is the mechanical energy, $M$ is the total mass in the system, $\phi_m$ is the mass flow, $\phi_w$ is a source of mechanical energy (i.e. work done on the system) and $\epsilon$ is the loss of mechanical energy (i.e. due to friction). $\Delta\phi_me$ is the difference in mechanical energy flowing in and out of the system.
The losses are generally modelled by the Darcy-Weisbach equation:
$$\epsilon=f_D\frac{L}{D}\frac{1}{2}v^2$$
where $f_D$ is the darcy friction factor, $L/D$ is the pipe length to diameter ratio and $v$ is the velocity. This should get you going on solving the problem with viscosity. However, i must stress that perhaps a CFD approach is more suited in this case than an analytical treatment.
Best Answer
There is a hypothesis left over here: the flow is irrotational, i.e. the so-called vorticity
$$\vec{\omega}=\nabla\times\vec{v}$$
is zero. As a result, the velocity field is a gradient,
$$\vec{v}=\nabla\phi,$$
for some scalar field $\phi$: the field appearing in the Bernoulli equation you quoted. Note that the right-hand side is not necessarily zero but it can be taken as a constant, which does not depend either on position or on time.
This is a simple consequence of (i) Euler equation for an incompressible fluid,
$$\frac{\partial\vec{v}}{\partial t} + (\vec{v}\cdot\nabla) v + \nabla\left(\frac{p}{\rho}+\psi\right)=0,$$
where $\psi=gh$ but it could be the potential for any other conservative force field; and (ii) the identity
$$\vec{v}\times\vec{\omega} = \nabla\left(\frac{v^2}{2}\right) - (\vec{v}\cdot\nabla) v.$$
Then Euler equation reads
$$\nabla\left(\frac{\partial\phi}{\partial t}+\frac{v^2}{2}+\frac{p}{\rho}+\psi\right)=0.$$
Therefore the term between parentheses is function of time only but we can absorb any time dependence into $\phi$, and therefore we end up with
$$\frac{\partial\phi}{\partial t}+\frac{v^2}{2}+\frac{p}{\rho}+\psi=C$$
where $C$ is constant over space and time. This is therefore a stronger result than the traditional Bernoulli theorem where the constancy is only along each streamline, with a priori a different constant for each streamline. Here the constancy is over the entire fluid.