Using TEBD2D for a long range 1D chain

[shadibeh]

Good time;
I would like to apply TEBD to compute two-point correlation function for a 1D chain (PBC) with three-body interactions. Motivated by @jcmgray's suggestion in #67, and quimb document of
https://quimb.readthedocs.io/en/latest/examples/ex_TEBD_evo.html)
I have tried to use TEBD2D through the following code:

import numpy as np
import quimb as qu
import quimb.tensor as qtn
L = 5 # number of qubits
zeros = '10' * (L//2)  #devide the vector in three pieces
binary = zeros   # the initial state
print('psi0:', f"|{binary}>")
# Initial state is turned into a matrix product state using MPS_computational_state():
psi0 = qtn.MPS_computational_state(binary)
psi0.show()  # prints ascii representation of state
I = qu.eye(2)
X = qu.pauli('X')
Z = qu.pauli('Z')
# the dense hamiltonian we are aiming for
H_dense = sum(
    qu.ikron([Z, X, Z], dims=[2] * L, inds=[i, i + 1, i + 2])
    for i in range(L - 2)
)

H_dense = H_dense+ sum(
    qu.ikron([Z, X, Z], dims=[2] * L, inds=[l-2, l-1, l-L])
    for l in range(L, L+1)
)

H_dense = H_dense + sum(
    qu.ikron([Z, X, Z], dims=[2] * L, inds=[k-1, k-L, k-L+1])
    for k in range(L, L+1)
)

H = qtn.LocalHam2D(L, 1, H_dense, H1=None)
# Setup TEBD object to find the evolved state
tebd = qtn.TEBD2D(psi0, H)
# Setup parameters
# times we are interested in
ts = np.linspace(0, 80, 101)
# compute the ZZ correlation on qubits 3 & 4
mzz_t_j = []
# range of bonds, and sites
js = np.arange(0, L)
bs = np.arange(1, L)
# generate the state at each time in ts
#     and target error 1e-3 for whole evolution
for psit in tebd.at_times(ts, tol=1e-3):
    mzz_j = []

    # there is one more site than bond, so start with mag
    #     this also sets the orthog center to 0
    mzz_j += [qu.pauli_correlation(psit, ss='ZZ', sysa=3, sysb=4)]
    mzz_t_j += [mz_j]

However, I get the following error: AttributeError: 'TEBD2D' object has no attribute 'at_times'
Would you please let me know how to resolve the issue?
Thanks a lot