def get_pij_method(model=F81, frequencies=None, kappa=None):
    """
    Returns a function for calculation of probability matrix of substitutions i->j over time t.

    :param kappa: kappa parameter for HKY model
    :type kappa: float
    :param frequencies: array of state frequencies \pi_i
    :type frequencies: numpy.array
    :param model: model of character evolution
    :type model: str
    :return: probability matrix
    :rtype: function
    """
    if is_f81_like(model):
        mu = get_mu(frequencies)
        return lambda t: get_f81_pij(t, frequencies, mu)
    if JTT == model:
        return get_jtt_pij
    if HKY == model:
        return lambda t: get_hky_pij(t, frequencies, kappa)