def METHOD_NAME(iels, rx, ree):
    rvals = eval_dofs(iels, rx)
    rng = rvals.max() - rvals.min()
    n_components = rvals.shape[-1]
    msd = 0.0
    for ic in range(n_components):
        rval = rvals[..., ic]
        sd = rval[:, ree]
        # ~ max. second derivative.
        msd += nm.abs(sd[..., 0] + sd[..., 2]
                      - 2.0 * sd[..., 1]).max(axis=-1)
    msd /= n_components
    return msd, rng