Here's my explanation for the algorithm for making the cards.
We can describe the cards as a matrix, with a row for each card and a column for each symbol. A cell will have a 1 if the card corresponding to the row has the symbol corresponding to the column, and 0 otherwise.
It turns out that for each prime number $p$ we can create such a solution, with $p^2+p+1$ cards and a total of $p^2+p+1$ symbols, with $p+1$ symbols on each card. The Dobble game has $p=7$, which gives 57 cards and 57 symbols with 8 symbols in each card. Here is the matrix for $p=5$, meaning there are 31 cards with 6 symbols in each:

Here's Python code for producing this solution:
def dobble(p):
cards = [[] for i in range(p**2 + p + 1)]
cards[0].append(0)
for i in range(p+1):
for j in range(p):
cards[1 + i*p + j].append(i)
cards[i].append(1 + i*p + j)
for i in range(p):
for j in range(p):
for k in range(p):
cards[1 + p + i*p + k].append(1 + p + j*p + (i*j - k) % p)
return cards
from itertools import combinations
for card0, card1 in combinations(dobble(7), 2):
assert len(set(card0) & set(card1)) == 1
The matrices $C_{ij}$ are all permutation matrices, that is, each row and each column have exactly one 1. They are all just a cyclically shifted reverse-diagonal, with the diagonal shifted by $ij \mod p$.
It is quite quite easy to see that any two rows have exactly one column with a common 1, except for two rows that come from different rows in the large matrix of matrices. Let's see why they have exactly one common column.
A matrix $C_{ij}$ has a 1 in position $k,l$ if $ij = k + l$. All the arithmetic is done modulu $p$. The first row can be described by $i_0$ and $k_0$ (k is the row in the $C_{i_0}j$ matrices), and the second row by $i_1$ and $k_1$. Since we come from different rows in the large matrix, $i_0 \ne i_1$. A column is described by $j$ and $l$. If both rows have a 1 in the column, it means that:
- $i_0 j = k_0 + l$
- $i_1 j = k_1 + l$
Subtracting the equations gives us:
$(i_1 - i_0)j = k_1 - k_0$
Now, since $p$ is prime, and $i_0 \ne i_1$, we can solve this equation, and get a single result for $j$, namely
$j = (k_1 - k_0)(i_1 - i_0)^{-1}$
which we can put back in equation 1 to get $l$. So there is exactly one column which has 1 in both rows.
So, thanks to the integers modulu $p$ being a field, we have exactly one symbol which appears in each pair of cards.
Notes:
- This would have also worked with regular diagonals, but this way the matrix is symmetric.
- I based my understanding on Paige L.J., Wexler Ch., A Canonical Form for Incidence Matrices of Projective Planes, which was referenced in this answer to a related question.
- The code for generating the figure is here.