Short summary: Your algorithm should be very efficient; it looks like it only has to explore tens of thousands of subcases, at most. See the second half of this answer for a detailed analysis. However, there are other more general approaches that can also be used for this problem, which I will list next. They may be of independent interest, because they can be used for other sorts of counting problems in this vein.
Have you tried using MINION, a solver intended for counting the number of solutions to a constraint satisfaction problem? It looks like your problem should be encodable in MINION's input language.
You could also try any #SAT solver. A #SAT solver counts the number of satisfying assignments to a particular formula $\varphi$, provided in SAT. It would be easy to write down a SAT formula on 100 unknowns that is satisfied exactly by the $10\times 10$ matrices with five zeros and five ones in each row and column. Then, you could try any off-the-shelf #SAT solver and see how efficient it is.
I don't know how well these will deal with the high amount of symmetry present in this problem, but it wouldn't be too hard to give them a try and see.
For more on these subjects, see the following two questions:
Another option, as suggested by Steven Stadnicki, is to express the set of valid matrices as a BDD. You should be able to construct a BDD on 100 variables that represents the set of all valid matrices, i.e., the set of satisfying assignments to the formula $\varphi$. Then, there are standard constructions to count the number of solutions to a BDD, using dynamic programming (see e.g., Knuth, vol. 4, TAOCP). This is a great suggestion. You'll have to try it to see how large a BDD you get; asymptotically, the size of the BDD will be exponential in the dimensions of the matrix, but for a $10 \times 10$ matrix, if you pick a good variable ordering, it's plausible that the BDD might be small enough that you can make this approach work for your particular problem. Credit to Steven Stadnicki for this suggestion.
Efficiency of your approach
You also asked about the efficiency of your approach, and in particular, the number of subcases it will explore. We can analyze this by looking at the number of partitions of an integer. The number of partitions of $n$ is denoted $p(n)$. For instance, $p(3)=3$, since the three partitions are $3=3$, $3=2+1$, $3=1+1+1$.
For starters, let's look at the number of subcases you examine for $M=6$ (i.e., when examining $6 \times 10$ matrices). With a bit of thinking, you can realize that each subcase corresponds to a unique partition of $30$. For instance, the subcase where you want column sums of $4,4,4,3,3,3,3,2,2,2$ (in your notation, $c_0=c_1=c_5=0$, $c_2=3$, $c_3=4$, $c_4=3$) corresponds to the partition $30=4+4+4+3+3+3+3+2+2+2$. Thus, the set of subcases is a subset of the set of partitions of $30$. (Not every partition of $30$ corresponds to a subcase; only those partitions where each summand in the partition is in the range 1..5.) Therefore, the number of subcases is $\le p(30)$. Consulting OEIS, we find that $p(30)=5604$.
It follows that for the case $M=6$, where you examine $6 \times 10$ matrices, your algorithm explores at most 5604 subcases.
This analysis can be generated to all $M \ge 6$. For $M=7$, it naively looks like we're looking at partitions of $35$: except that each summand must be in the range 2..5, and there must be exactly 10 summands. If we subtract one from each summand, we get a partition of $25$ (where each resulting summand is in the range 1..4). Therefore, the number of subcases for $M=7$ is $\le p(25)=1958$. By the same reasoning, the number of subcases for $M=8,9,10$ is at most $p(20),p(15),p(10)$, i.e., at most $627, 176, 42$, respectively.
Next, let's look at $M=5$. Now the column sums are in the range 0..5, and they must sum to $25$. One loose bound can be obtained by adding one to each column sum, and noting that each subcase corresponds to a partition of $35$ where all summands are in the range 1..6; there are thus $\le p(35)=14883$ such subcases. A slightly better bound can be obtained by counting the number $c_0$ of columns whose column sums are $0$; for each fixed $c_0$, the number of subcases is $\le p(25)$ (since the other column sums sum to $25$). The number $c_0$ can be anything in the range 0..5. Therefore, the total number of subcases is $\le 6 \times p(25) = 11748$.
What about $M=4$? Here the column sums are each in the range 0..5, and they must sum to $20$. Repeating the above reasoning gives one upper bound $\le p(30) = 5604$. The slightly better upper bound is $\le 7 \times p(20) = 4389$.
Continuing in this way, for $M=3,2,1,0$ the number of subcases is at most $1408,378,70,1$. I'm sure you can get better bounds for many of these cases, but this is sufficient to give an upper bound on the number of subcases.
Summing this up over all $M$, we find that the total number of subcases is at most
$$1+70+387+1408+4389+11748+5604+1958+627+176+42 = 26410.$$
Therefore, your algorithm explores at most 26410 subcases. (This is a conservative upper bound and the true number will be lower than this.) This is small enough that your algorithm should be very efficient, for the case of $10 \times 10$ algorithms.