Based on your description of the problem, I'm going to assume that you're after the first moment of a *Skellam distribution truncated at zero* and not the partial first moment of a Skellam (which is the expression you give in the title of your question). **If my assumption is wrong** and you do want a partial moment, the answer is similar with a slight modification, see the end of this answer before the code block.

Then, changing notation slightly, the pmf is:
$$
p(u; \lambda_1, \lambda_2) = (\lambda_1/\lambda_2)^{u/2} I_u(2\sqrt{\lambda_1 \lambda_2}) Z_{\lambda_1, \lambda_2}^{-1}
$$

Where $Z_{\lambda_1, \lambda_2} = \sum_{k=0}^\infty (\lambda_1/\lambda_2)^{k/2} I_k (2\sqrt{\lambda_1\lambda_2})$ is the normalizing constant, after cancelling out $e^{-\lambda_1-\lambda_2}$.

I follow the same approach as in my answer to this question on Cross-Validated. Briefly,

The moment generating function (mgf) is,
$$
\begin{aligned}
\mathcal{M}(t; \lambda_1, \lambda_2) &= Z_{\lambda_1, \lambda_2}^{-1} \sum_{u=0}^\infty e^{tu} (\lambda_1/\lambda_2)^{u/2} I_u (2\sqrt{\lambda_1\lambda_2}) \\
&= Z_{\lambda_1, \lambda_2}^{-1} \sum_{u=0}^\infty (\frac{\lambda_1 e^t}{\lambda_2 e^{-t}})^{u/2} I_u (2\sqrt{\lambda_1 e^t \lambda_2 e^{-t}}) \\
&= Z_{\lambda_1, \lambda_2}^{-1} Z_{\lambda_1 e^t, \lambda_2 e^{-t}}
\end{aligned}
$$

In practice the normalizing constant can be calculated from the easy-to-compute special function Marcum's Q (or, equivalently, from the noncentral-$\chi^2$ cdf, see code block below),
$$
Z_{a,b} = Q(\sqrt{2b},\sqrt{2a}) e^{a+b}
$$

Differentiating the mgf around $t=0$ and simplifying gives the expectation of the truncated distribution:
$$
\begin{aligned}
\mathcal{M}'(t; \lambda_1, \lambda_2)|_{t=0} &= \mathbb{E}[u; \lambda_1, \lambda_2] \\
&= (\lambda_1 - \lambda_2) + Z_{\lambda_1,\lambda_2}^{-1} (\lambda_2 I_0 (2\sqrt{\lambda_1 \lambda_2}) + \sqrt{\lambda_1 \lambda_2} I_1(2\sqrt{\lambda_1\lambda_2}))
\end{aligned}
$$

**If you are actually after the partial first moment**, just multiply the expectation above through with the normalizing constant (and cancelled terms):
$$
\mathbb{E}_{u \geq 0}[u; \lambda_1, \lambda_2] = \mathbb{E}[u; \lambda_1, \lambda_2] \cdot Z_{\lambda_1,\lambda_2} e^{-\lambda_1 - \lambda_2}
$$

Example in R:

```
Q <- function(b,a) 1-pchisq(b^2, 2, a^2)
Z <- function(a,b) Q(sqrt(2*b),sqrt(2*a)) * exp(a+b)
mean_u <- function(l1,l2)
(l1 - l2) + (l2*besselI(2*sqrt(l1*l2),0) +
sqrt(l1*l2)*besselI(2*sqrt(l1*l2),1))/Z(l1,l2)
l1 <- 4; l2 <- 12
mean_u(l1,l2)
# [1] 0.9668963
# compare to naively computing series,
u_grid <- 0:100
pmf_u <- sapply(u_grid, function(u)
exp(-l1-l2)*(l1/l2)^(u/2) * besselI(2*sqrt(l1*l2),u))
sum(u_grid * pmf_u)/sum(pmf_u)
# [1] 0.9668963
# or (inefficiently) simulating via rejection sampling
set.seed(1)
k <- rpois(10^7,l1) - rpois(10^7,l2)
mean(k[k>=0])
# [1] 0.9637044
## alternatively, partial first moment
mean_u(l1,l2) * Z(l1,l2) * exp(-l1-l2)
# [1] 0.02553668
sum(u_grid * pmf_u) #naive summation
# [1] 0.02553668
```