T = 100000; t = 0; lambda = 2; mu = 0.25; X = 0; jump_times = [0]; jump_chain = [X]; while t <= T t = t - log(rand) / (lambda + X*mu); if t > T break; else jump_times = [jump_times t]; if rand < lambda / (lambda + X*mu) X = X + 1; else X = X - 1; end jump_chain = [jump_chain X]; end end