T = 10; t = 0;

Q = [-9 2 3 4; 1 -4 3 0; 1 2 -7 4; 1 0 3 -4];
J = [0 2/9 3/9 4/9; 1/4 0 3/4 0; 1/7 2/7 0 4/7; 1/4 0 3/4 0];

U = rand; X = (U <= 1/2) + 2*(U > 1/2);
jump_times = [0]; jump_chain = [X];

while t <= T
    t = t + log(rand) / Q(X,X);
    if t > T
        break;
    else
        jump_times = [jump_times t];
        X = min(find(rand<cumsum(J(X,:))));
        jump_chain = [jump_chain X];
    end
end