monte carlo - How to simulate a jump-diffusion process?

I would like to price Asian and Digital options under Merton's jump-diffusion model. To that end, I will have to simulate from a jump diffusion process.

In general, the stock price process is given by

$$S(t) = S(0)e^{(r-q-\omega)t+X(t)},$$ where $\omega$ is the drift term that makes the discounted stock price process a martingale and where $X(t) = \sigma W(t)+\sum_{i=1}^{N(t)} Y_i$ is a jump-diffusion process. The process $N(t)$ is a Poisson process with intensity $\lambda$ and the jump sizes $Y_i$ are iid Normally distributed $Y_i \sim N(\mu, \delta^2)$.

So, the issue here is the simulation of $X(t)$. I would like to simulate $X(t)$ as follows: first, I construct a grid $[0,T]$ with grid step size $\Delta t = 1/252$ (number of trading days in one year). I now want to compute $X$ on this grid stepwise. Therefore:

$$X_{n\Delta t} = X_{(n-1) \Delta t}+\sigma \sqrt{\Delta t} \epsilon_n + p_n,$$ where $\epsilon_n$ are standard normal random variables (easy to simulate in Matlab) and where $p_n$ are Compound Poisson random variables independent of $\epsilon_n$. The question is, how to simulate from these $p_n$?

I was thinking to write a compound Poisson generator myself:

function [p] = CPoissonGenerator(mu,delta,lambda,dt,rows,columns)
   N = poissrnd((1/lambda)*dt,rows,columns);
   Y = normrnd(mu,delta,1,N);
   cP = cumsum([0,Y]);
   p = cP(end);
end

However, I do not know how to generate a matrix of Poisson random variables without using a loop in the above code? This can be very time consuming. Furthermore, I'm not sure if the program is correct. Is there a better way to do this?