Converting logit coefficients into probabilities through the inverse multinomial logit

I am trying to come up with the transition intensities of a parameterized Markov model, i.e, the elements of the transition matrix weighted by the covariates through the inverse Logit relation:

$$\Theta= \begin{bmatrix} \theta_{11}&\theta_{12}\\ \theta_{21}&\theta_{22} \end{bmatrix}$$

where as usual, the matrix sums row-wise should equal 1, and the $\theta_{ij}$ are calculated thus (let's say we have as explanatory variables an intercept, Income, and a Location dummy taking values 0 or 1):

$$\theta_{ij}=\frac{e^{\beta^{ij}_{0}+\beta^{ij}_{1}Income+\beta^{ij}_{2}LocationDummy}}{1+\sum_\limits{i\neq j}e^{\beta^{ij}_{0}+\beta^{ij}_{1}Income+\beta^{ij}_{2}LocationDummy}}$$

I've already estimated the relevant parameters through Maximum Likelihood and now need to find the $\theta_{ij}$, denoting the transition probability between states $i$ and $j$ as per the formula above to fill in the transition matrix.

The trouble is that for most values of $\beta$'s times the covariates that I include, the expression goes to either $0$ or $+\infty$. If, say, income is 500, then the expression can very easily go to one of these bounds, whereas it should be a number between 0 and 1.

I have estimated the betas, and am now plugging them in the expression for $\theta_{ij}$ to get the matrix $\Theta$.

If it helps, I am using R, and know the invlogit function, which does essentially this, but for the standard logit. Is there some function in R that does this, or does this issue have some canonical fix, like normalizing income values?