We consider the random Hermitian matrix model of dimension $2n$, with external source, defined by the probability density function \begin{equation*} \frac{1}{Z_{2n}} \lvert \det(M) \rvert^α e^{-2n\mathrm{Tr} (V(M) - AM)}, \quad V(x) = \frac{x^4}{4} - t\frac{x^2}{2}, \end{equation*} where the external source $A$ has two eigenvalues $\pm a$ of equal multiplicity. We investigate the limiting local statistics of the eigenvalues of $M$ around $0$ in certain critical regimes as $n \to \infty$. When the parameters $t$ and $a$ lie on a critical curve along which the limiting mean eigenvalue density vanishes as $|x|^{1/3}$, the double scaling limit of the correlation kernel is constructed from functions associated with the Boussinesq equation. This new limiting kernel reduces to the classical Pearcey kernel when $α= 0$. Furthermore, in the multi-critical case where the limiting mean eigenvalue density vanishes as $|x|^{5/3}$, the limiting kernel is built from the second member of the Boussinesq hierarchy. We derive the results by transforming the random matrix model into biorthogonal ensembles that are analogous to the Muttalib-Borodin ensemble, and then analyzing its asymptotic behavior via a vector Riemann-Hilbert problem.