The delta method states
$$ \operatorname{Var}(g(X)) = [g'(X)]^2 \operatorname{Var}(X)$$
Because this problem involves two parameters, we can extend this to the multivariate delta rule
$$ =\nabla g(X)^T \, \Sigma \, \nabla g $$
Here,
$$ g = \left[e^{\beta_{0}+\beta_{1}}-e^{\beta_{0}}\right] /\left[\left(1+e^{\beta_{0}+\beta_{1}}\right)\left(1+e^{\beta_{0}}\right)\right] $$
and $\Sigma$ is the variance covariance matrix from your model. $\nabla g$ is...gross. I'm not going to do that by hand, and computer algebra while fast yields a mess of symbols. You can however use autodifferentiation compute the gradient. Once you calculate the variance, then its simply your estimate of the difference in probs plus/minus 1.96 times the standard deviation (root of the variance). Caution, this approach will yield answers below 0 or above 1.