%% computes the transition matrix of the model given parameters and current state
function eQ = exQ(nr, nl, N, params)
    arl = params(1);
    brl = params(2);
    crl = params(3);
    ars = params(4);
    brs = params(5);
    crs = params(6);
    asr = params(7);
    bsr = params(8);
    csr = params(9);
    Q = zeros([3 3 length(nr)]);
    Q(1,2) = arl + brl.*nl/N + crl.*(N-nr-nl)/N;
    Q(1,3) = ars + brs.*(N-nr-nl)/N + crs.*nl/N;
    Q(1,1) = -Q(1,2)-Q(1,3);
    Q(2,1) = arl + brl.*nr/N + crl.*(N-nr-nl)/N;
    Q(2,3) = ars + brs.*(N-nr-nl)/N + crs.*nr/N;
    Q(2,2) = -Q(2,1)-Q(2,3);
    Q(3,1) = asr + bsr.*nr/N + csr.*nl/N;
    Q(3,2) = asr + bsr.*nl/N + csr.*nr/N;
    Q(3,3) = -Q(3,1)-Q(3,2);
    eQ = expm(Q);
end