Abstract indices is a good alternative to work with, just don't be afraid of using different ranges of symbols dedicated to each bundle you deal with. For instance, I prefer to use the letters $a$, $b$, $c$ etc (the initial segment of the Latin alphabet) as the labels for the tangent and cotangent bundles (depending on the position), and juxtapositions of these letters denote various tensor bundles. Using parentheses and brackets we can indicate tensor parts of these bundles.
Thus, as you correctly noticed $E_{a[bc]}$ can be seen as the bundle of 2-form-valued 1-forms, and this notation is indeed standard. If there is another bundle that I want to be treated differently from the tensor bundles, I would use letters $\Phi$, $\Psi$, $\Xi$ etc associated to this bundle, and again tensor products of thus bundle and its dual will be simply denoted by using letter $E$ with the letters from the specified range attached as indices in various positions. We may have something like $E^{\Phi \Psi}{}_{\Xi}$ or $E^{[\Phi \Psi]}{}_{\Xi}$, and so on. The question is how to introduce the usual operations in this notation. This has been already done. See the references below.
For instance, the exterior covariant derivative of a 1-form with values in 2-forms will be written as simple as $2\,\nabla_{[i}P_{j][bc]}$. The notation $P_{[a}Q_{b]}$ can be used for the wedge product $P\wedge Q$ (up to the coefficient). A connection in a vector bundle $E^{\Phi}$ will be then a map $\nabla_{a}\colon E^{\Phi}\to E^{\Phi}{}_{a}$ admitting a certain sloppiness of the language. Another useful trick with the abstract indices is not always to use the indices ("suppressing") which comes naturally when one calculates in this notation a lot. For the Hodge dual you may probably want to use the volume form $\varepsilon_{a^{1}\dots a^{n}}$ which is a $n$-form with certain properties.
References:
- R.Penrose, W.Rindler "Spinors and space-time. Vol.1" which is usually referred to as the origin of the abstract index notation (in fact, it was known to physicists a while before the book appeared, but R.Penrose gave a logical foundation to expose its mathematical rigor in full glory).
- R.Wald "General Relativity", see Appendix B where one can find useful formulae for differential forms given in the abstract index notation. This book is somewhat more accessible then the previous one.
- R.W.R.Darling "Differential Forms and Connections" is more suitable for mathematicians, and is a useful reading as well. The author uses the classical "invariant" notation, such as $d^{\nabla}\omega$ for exterior covariant derivatives.
Edit. As requested I attempt to apply the abstract index notation to show the Leibniz rule for the case of a 1-form $v_{a}$ and a 2-form $w_{ab} = w_{[ab]}$. So, we want to show
$$
(d(v \wedge w))_{abcd} = (dv \wedge w)_{abcd} + (-1)^{deg(v)}(v \wedge dw)_{abcd}
$$
The definitions disentangled give
$$
(v \wedge w)_{bcd} = \frac{3!}{1! \cdot 2!} v_{[b} w_{cd]}
$$
$$
(d(v \wedge w))_{abcd} = 4 \nabla_{[a} (v \wedge w)_{bcd]} = \frac{4!}{1! \cdot 2!} \nabla_{[a} (v_{b} w_{cd]}) \tag{1}
$$
$$
dv_{ab}=2 \nabla_{[a}v_{b]}
$$
$$
(dv \wedge w)_{abcd} = \frac{4!}{2!\cdot 2!} (2 \nabla_{[a}v_{b})w_{cd]} = \frac{4!}{1!\cdot 2!} (\nabla_{[a}v_{b})w_{cd]} \tag{2}
$$
$$
(dw)_{bcd}=3 \nabla_{[b} w_{cd]}
$$
$$
(v \wedge dw)_{abcd} = \frac{4!}{1! \cdot 3!} (3 v_{[a} \nabla_{b} w_{cd]}) = \frac{4!}{1! \cdot 2!} ( v_{[a} \nabla_{b} w_{cd]}) \tag{3}
$$
It remains to compare (1) with (2) and (3) to understand what is going on: the alternation is responsible for the sign $(-1)^{deg(v)}$ which is just $-1$ in our case.
Remark. If our forms have coefficients in vector bundles we need to take care to define what the wedge product really means for the calculations to make sense.
I think the answers to your first two questions are yes, there is always a codifferential; and no, there is not always a Hodge decomposition. I don't know anything about the Einstein manifold case.
Let's say $E$ is a vector bundle over $M$ with metric and compatible connection $\nabla$. ($E$ could be some tensor bundle with $\nabla$ induced from the Levi-Civita connection, for example.) As you say, $\nabla$ gives an exterior derivative $d^\nabla$ on $E$-valued differential forms that obeys the Leibniz rule
$$ d^\nabla (\omega \otimes e) = d\omega \otimes e + (-1)^{\text{deg} \omega} \omega \wedge \nabla e ,$$
where $\omega$ is a homogeneous form and $e$ is a section of $E$, and $\omega \wedge \nabla e$ means $\sum_i (\omega \wedge dx^i) \otimes \nabla_{\partial_i} e$.
For each degree of forms $k$, you can view $d^\nabla$ as a map $d^\nabla: C^\infty( \Lambda^k \otimes E) \to C^\infty( \Lambda^{k+1} \otimes E)$. Using the metrics, each of those spaces of sections are endowed with $L^2$ inner products. Then $d^\nabla$ always has a formal $L^2$-adjoint $\delta^\nabla$ going the other way: $\delta^\nabla: C^\infty( \Lambda^{k+1} \otimes E) \to C^\infty( \Lambda^k \otimes E)$.
(Note the space of smooth sections can be completed to the space of $L^2$ sections, but $d^\nabla$ and $\delta^\nabla$ are unbounded operators and can generally only be defined on some dense subspace of that $L^2$ space.)
So you can always define $d^\nabla$ and its formal adjoint $\delta^\nabla$. The issue is that $d^\nabla$ squares to zero if and only if the connection $\nabla$ on $E$ is flat, meaning its curvature is zero. (Flat vector bundles $E$, i.e., bundles on which there exists a flat connection, are relatively scarce.) We need $d^\nabla$ to square to zero to even attempt to define the de Rham cohomology of $E$-valued forms, and to have nice properties like the Hodge decomposition. (I want to say that the "Dirac operator" $D = d^\nabla + \delta^\nabla$ and the "Laplacian" $\Delta = D^2$ are not elliptic unless $d^\nabla$ squares to zero, but I need to think about that.)
I don't have a good reference for this at hand, although the Wikipedia article on bundle-valued forms might be useful.
Best Answer
I am aware of (and tried to work with) the following three options:
Personally, I find Cadabra the most suitable for my current needs (extensive calculations with polynomial tensor expressions), but they differ from yours. Cadabra is quite good at taking tensor products and dealing with the product rule when differentiating. It can distinguish ranges of indices to represent different bundles. Lots of things can be done by making substitution rules. It also has a lot of machinery that I have never used (e.g., for tensor symmetries). It can work with metrics nicely, by the way.
I'm curious too about other possible pieces of software, other people experience, and success stories.