Since there are no free extraneous charges anywhere
$$div \vec {D} = \dfrac {\partial D_{x}}{\partial x} = 0$$ or, $$D_{x} = Constant$$
But $$D_{x} = 0$$ at $$\infty$$, so, $$D_{x} = 0$$, every where.
Thus, $$\vec {E} = -\dfrac {\vec {P}_{0}}{\epsilon_{0}} \left (1 - \dfrac {x^{2}}{d^{2}}\right )$$ or, $$E_{x} = -\dfrac {P_{0}}{\epsilon_{0}} \left (1 - \dfrac {x^{2}}{d^{2}}\right )$$
So, $$\varphi = \dfrac {P_{0}x}{\epsilon_{0}} - \dfrac {P_{0}x^{3}}{3\epsilon_{0}d^{2}} + constant$$
Hence,
$$\varphi (+d) - \varphi (-d) = \dfrac {2P_{0}d}{\epsilon_{0}} - \dfrac {2P_{0}d^{3}}{3d^{2}\epsilon_{0}} = \dfrac {4P_{0}d}{3\epsilon_{0}}$$.