It should be possible to restate that as $P(\mu-\sigma \Phi^{-1}(\frac{p+1}{2})\leq X\leq \mu+\sigma \Phi^{-1}(\frac{p+1}{2}))=p$.
In this answer, it says:
For a normal distribution, the probability of being within $\Phi^{-1}\left(\frac{p +1}{2}\right)$ standard deviations of the mean is $p$, where $\Phi^{-1}$ is the inverse of the cumulative distribution of a standard normal.
I tried expressing $\Phi^{-1}\left(\frac{p +1}{2}\right)$ in terms of $\operatorname{erf^{-1}}$, but then again I can't get rid of the error function.
Also taking $\Phi$ on both sides would give $\Phi(p)=(p+1)/2$, but a simulation with MATLAB for the case described in the linked question shows it checks out.
(All this provided I interpreted the linked answer correctly.)
If $X\sim N(\mu,\sigma)$, then $Y=\frac{X-\mu}{\sigma}\sim N(0,1) $ and
$$ \mathbb{P}[\mu-k\sigma \leq X \leq \mu+k\sigma] = \mathbb{P}[-k\leq Y \leq k].$$ Can you recognize $\Phi$ in the RHS now?