def Jwik(dW, h, n=5):
    """matrix J approximating repeated Stratonovich integrals for each of N
    time intervals, using the method of Wiktorsson (2001).

    Args:
      dW (array of shape (N, m)): giving m independent Weiner increments for
        each time step N. (You can make this array using sdeint.deltaW())
      h (float): the time step size
      n (int, optional): how many terms to take in the series expansion

    Returns:
      (Atilde, J) where
        Atilde: array of shape (N,m(m-1)//2,1) giving the area integrals used.
        J: array of shape (N, m, m) giving an m x m matrix of repeated
        Stratonovich integral values for each of the N time intervals.
    """
    m = dW.shape[1]
    Atilde, I = Iwik(dW, h, n)
    J = I + 0.5*h*np.eye(m).reshape((1, m, m))
    return (Atilde, J)