It is perhaps not so surprising that modular forms mod p do not always behave the way one would expect based on properties of modular forms in characteristic 0. Here is an example of such a discrepancy. (This is of course "well-known", and as with a lot of "well-known" things, it's not always easy to find a reference.)

It is a classical result that the space of modular forms M:=M_k(\textrm{SL}_2(\mathbb{Z})) has a basis of Hecke eigenforms. (See for instance Section 4.2.3 in Kilford's Modular forms: a classical and computational introduction.) If you combine this with the q-expansion principle, you get a multiplicity one result, i.e. the eigenspace associated to any system of Hecke eigenvalues occurring in M is one-dimensional. There are similar results for level N, where one needs to differentiate between oldforms and newforms. (See for instance Section 5.8 in Diamond-Shurman's A first course in modular forms.)

Let's now see what happens when we instead work with \tilde{M}:=\tilde{M}_k(\textrm{SL}_2(\mathbb{Z})), the space of modular forms mod p. It turns out that this space need not have a basis of Hecke eigenforms. Here is a sketch of proof: suppose these spaces have a basis of Hecke eigenforms, then each eigensystem has multiplicity one by the q-expansion principle. But a theorem of Serre-Tate says that there are only finitely many Hecke eigensystems mod p (all weights put together). Say this number is H. Choosing a weight k such that the space of modular forms of weight k has dimension greater than H gives us more than H eigensystems in weight k, contradiction.

The argument tells us where to look for such situations. I'll restrict my attention to cusp forms, but this is not important. If p is 5, then there are only 4 Hecke eigensystems in all weights. So taking a space of dimension at least 5 would do the trick, and indeed the space of cusp forms of weight 60 works. But it is possible to find smaller such examples (of course, the space must be at least two-dimensional for this to work). For p=5, the earliest example comes up in weight 36. Sage makes it easy to explore this as follows. First, we need a basis for the space of cusp forms mod 5 of weight 36, and we'll use the corresponding Victor Miller basis:

sage: bas = victor_miller_basis(k=36, prec=20, cusp_only=True)
sage: bas_mod5 = [ f.change_ring(GF(5)) for f in bas ]

The second line reduces the original Victor Miller basis (which lives over the integers) modulo 5. Now we figure out the matrix of the Hecke operator T_2 on this basis. This should just be a matter of doing hecke_operator_on_basis(bas_mod5, 2, 36), but at the moment that doesn't work over finite fields (for silly reasons, and this should be an easy bug to fix).

No matter, we can just do it by hand:

sage: f1, f2, f3 = bas_mod5; print f1; print f2; print f3
q + 3*q^4 + 2*q^6 + 2*q^9 + 2*q^11 + O(q^12)
q^2 + 4*q^4 + 2*q^5 + q^6 + 4*q^7 + 4*q^8 + 3*q^9 + 2*q^10 + O(q^12)
q^3 + 3*q^4 + 4*q^5 + 2*q^6 + q^7 + 3*q^8 + q^9 + 4*q^10 + O(q^12)
sage: hecke_operator_on_qexp(f1, 2, 36)
q^2 + 2*q^3 + O(q^6)

And herein lies the beauty of the Victor Miller basis: we immediately see that the last result is really f_2+2f_3. Similarly:

sage: hecke_operator_on_qexp(f2, 2, 36)
q + 4*q^2 + q^3 + 2*q^4 + 2*q^5 + O(q^6)
sage: hecke_operator_on_qexp(f3, 2, 36)
3*q^2 + 2*q^3 + 3*q^4 + 4*q^5 + O(q^6)

Now we can write down the matrix of T_2 and see what happens:

sage: t2 = matrix(GF(5), [[0, 1, 0], [1, 4, 3], [2, 1, 2]]); t2
[0 1 0]
[1 4 3]
[2 1 2]
sage: t2.jordan_form()
[4|0 0]
[-+---]
[0|1 1]
[0|0 1]

And there you have it: since the matrix is not diagonalisable, this space of cusp forms cannot even have a basis of eigenvectors for T_2, let alone a basis of eigenvectors for all the Hecke operators.

Another "small" example is p=67, k=32, where the space of cusp forms is two-dimensional and the matrix of T_2 has Jordan normal form

sage: t2 = matrix(GF(67), [[0, 1], [5, 28]])
sage: t2.jordan_form()
[14  1]
[ 0 14]