<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def eigh(a, eigvec=True, rcond=None):
 """ Eigenvalues and eigenvectors of symmetric matrix ``a``. Args: a: Two-dimensional, square Hermitian matrix/array of numbers and/or :class:`gvar.GVar`\s. Array elements must be real-valued if `gvar.GVar`\s are involved (i.e., symmetric matrix). eigvec (bool):
 If ``True`` (default), method returns a tuple of arrays ``(val, vec)`` where ``val[i]`` are the eigenvalues of ``a`` (in ascending order), and ``vec[:, i]`` are the corresponding eigenvectors of ``a``. Only ``val`` is returned if ``eigvec=False``. rcond (float):
 Eigenvalues whose difference is smaller than ``rcond`` times their sum are assumed to be degenerate (and ignored) when computing variances for the eigvectors. Default (``rcond=None``) is ``max(M,N)`` times machine precision. Returns: Tuple ``(val,vec)`` of eigenvalues and eigenvectors of matrix ``a`` if parameter ``eigvec==True`` (default). The eigenvalues ``val[i]`` are in ascending order and ``vec[:, i]`` are the corresponding eigenvalues. Only the eigenvalues ``val`` are returned if ``eigvec=False``. Raises: ValueError: If matrix is not square and two-dimensional. """

    a = numpy.asarray(a)
    if a.dtype != object:
        val, vec = numpy.linalg.eigh(a)
        return (val, vec) if eigvec else val
    amean = gvar.mean(a)
    if amean.ndim != 2 or amean.shape[0] != amean.shape[1]:
        raise ValueError('bad matrix shape: ' + str(a.shape))
    if rcond is None:
        rcond = numpy.finfo(float).eps * max(a.shape)
    da = a - amean
    val0, vec0 = numpy.linalg.eigh(amean)
    val = val0 + [
        vec0[:, i].conjugate().dot(da.dot(vec0[:, i])) for i in range(vec0.shape[1])
        ]
    if eigvec == True:
        if vec0.dtype == complex:
            raise ValueError('cannot evaluate eigenvectors when a is complex')
        vec = numpy.array(vec0, dtype=object)
        for i in range(len(val)):
            for j in range(len(val)):
                dval = val0[i] - val0[j]
                if abs(dval) < rcond * abs(val0[j] + val0[i]) or dval == 0.0:
                    continue
                vec[:, i] += vec0[:, j] * (
                    vec0[:, j].dot(da.dot(vec0[:, i])) / dval
                    )
        return val, vec
    else:
        return val