There's another way to perform calculations in Rijndael's Galois field $GF(2^8)$ as shown in this video by Creel. It's less optimized in terms of computational effort but stays closer to the math mentioned in the Wikipedia article.
The following method works by transforming the bytes into polynomials, then multiplying them, and converting the resulting polynomial back to a byte. It's easier than it sounds as those conversions are pretty mechanical.
Let's use the numbers from your question (which also are found in the NIST publication on AES and this visualization of Rijndael at formaestudio.com):
$$
\begin{bmatrix}
\mathtt{02} & \mathtt{03} & \mathtt{01} & \mathtt{01} \\
\mathtt{01} & \mathtt{02} & \mathtt{03} & \mathtt{01} \\
\mathtt{01} & \mathtt{01} & \mathtt{02} & \mathtt{03} \\
\mathtt{03} & \mathtt{01} & \mathtt{01} & \mathtt{02} \\
\end{bmatrix}
\cdot
\begin{bmatrix}
\mathtt{d4} \\
\mathtt{bf} \\
\mathtt{5d} \\
\mathtt{30} \\
\end{bmatrix}
=
\begin{bmatrix}
\mathtt{04} \\
\mathtt{66} \\
\mathtt{81} \\
\mathtt{e5} \\
\end{bmatrix}
$$
The matrix has to go to the left of the column so that matrix multiplication is defined.
In this answer, I'll show the steps for multiplying the first column of the state with the first row of the fixed MixColumns matrix:
$$(\mathtt{d4} \cdot \mathtt{02}) + (\mathtt{bf} \cdot \mathtt{03}) + (\mathtt{5d} \cdot \mathtt{01}) + (\mathtt{30} \cdot \mathtt{01}) = \mathtt{04}$$
The centered dot $\cdot$ means multiplication over the finite field $GF(2^8)$, which ensures the results of our multiplications will fit into eight bits, one byte.
Addition corresponds to XOR in $GF(2^8)$.
First step: $\mathtt{d4}_{16} \cdot \mathtt{02}_{16}$
Initially, convert the hex representations of our bytes to binary. You can do this with bc but keep in mind that bc uses upper-case letters for hex digits:
bc <<< "obase=2;ibase=16; D4"
This gives us $\mathtt{11010100}_2$ for $\mathtt{d4}_{16}$.
$\mathtt{2}_{16}$ is $\mathtt{10}_2$.
We'll convert those two binary numbers to polynomials and multiply them, modulo 2.
For $\mathtt{11010100}_2$ we get this polynomial:
$\mathtt{x^7 + x^6 + x^4 + x^2}$
How to convert bits to polynomials? You take the indices of the binary number where the bit is 1 and use those indices as the exponents in the polynomial.
Notice that the exponents in the polynomial correspond to the indices in the binary number where the bit is set to 1:
$
\begin{array}{c|ccccccc}
bit & 1 & 1 & 0 & 1 & 0 & 1 & 0 & 0 \\
\hline
index & \color{red}{7} & \color{red}{6} & 5 & \color{red}{4} & 3 & \color{red}{2} & 1 & 0
\end{array}
$
In $\mathtt{11010100}$, the bit at index zero is 0, so exponent 0 doesn't go in the polynomial. Same for index one, its bit is 0. But bit at index two is 1, same for the bits at indices 4, 6, and 7. Thus, we create a polynomial with the exponents 7, 6, 4, and 2: $\mathtt{x^\color{red}{7} + x^\color{red}{6} + x^\color{red}{4} + x^\color{red}{2}}$
The polynomial for $\mathtt{10}_2$ is $\mathtt{x^1}$ since the bit at index one is 1.
Alright, we converted
$\mathtt{d4}_{16} \cdot \mathtt{02}_{16}$
to
$\mathtt{(x^7 + x^6 + x^4 + x^2) * x^1}$
Now we solve this polynomial multiplication modulo 2.
One way to calculate the product of polynomials is to write out the exponents of the two polynomials in a grid:
$
\begin{array}{c|c}
+ & 1 \\
\hline
2 & \\
4 & \\
6 & \\
7 & \\
\end{array}
$
And then add the exponents:
$
\begin{array}{c|c}
+ & 1 \\
\hline
2 & 3\\
4 & 5\\
6 & 7\\
7 & 8\\
\end{array}
$
Because we are in modulo 2, we need to remove all those exponents that occur an even number of times (two times, four times, etc.). But since the resulting exponents 3, 5, 7, and 8 only occur once—and 1 is an odd number—that's not necessary with this particular result.
We use those exponents to create the polynomial: $\mathtt{x^8 + x^7 + x^5 + x^3}$.
Now, convert the polynomial to a binary number: If there's an exponent in the polynomial, it means the bit at that index is 1:
$\mathtt{110101000}_2$
You can see that the bits at indices 8, 7, 5, and 3 are 1. Those indices correspond to the exponents of the polynomial $\mathtt{x^8 + x^7 + x^5 + x^3}$.
In hex, that binary number is $\mathtt{1a8}_{16}$ which is more than a byte, meaning the end result (after doing steps 2 to 5), too, might be larger than a byte which wouldn't work. In MixColumns, we transform one byte into another byte, not one byte into two bytes.
For these cases, the creators of Rijndael defined the polynomial $\mathtt{x^8+x^4+x^3+x^1+x^0}$ whose purpose is to reduce those intermediate results to a number equal to or smaller than a byte. This reducing polynomial is called $\mathtt{m(x)}$ in the literature. In other words, we subtract $\mathtt{m(x)}$ from our intermediate result $\mathtt{x^8 + x^7 + x^5 + x^3}$ to ensure the final multiplication result stays inside the finite field $GF(2^8)$.
Converting the reducing polynomial $\mathtt{x^8+x^4+x^3+x^1+x^0}$ to bits, we get $\mathtt{100011011}_2$. That's $\mathtt{11b}_{16}$ as, for example, calculated using printf and Bash' arithmetic expansion:
printf '%x\n' $((2#100011011))
Since subtraction corresponds to XOR in $GF(2^8)$, we XOR $\mathtt{1a8}_{16}$ with $\mathtt{11b}_{16}$:
printf "%x\n" $((0x1a8 ^ 0x11b))
Which yields a single byte: $\mathtt{b3}_{16}$. That's the result of $\mathtt{d4}_{16} \cdot \mathtt{02}_{16}$ in Rijndael's Galois field $GF(2^8)$.
Second step: $\mathtt{bf}_{16} \cdot \mathtt{03}_{16}$
Multiplying by $\mathtt{03}_{16}$ works the same as multiplying by $\mathtt{02}_{16}$:
Convert $\mathtt{bf}_{16} \cdot \mathtt{03}_{16}$ to binary:
$\mathtt{10111111}_2 \cdot \mathtt{11}_2$
Convert that to polynomials:
$\mathtt{(x^7 + x^5 + x^4 + x^3 + x^2 + x^1 + x^0) * (x^1 + x^0)}$
Solve the polynomial multiplication modulo 2:
$
\begin{array}{c|cc}
+ & 0 & 1 \\
\hline
0 & 0 & 1 \\
1 & 1 & 2 \\
2 & 2 & 3 \\
3 & 3 & 4 \\
4 & 4 & 5 \\
5 & 5 & 6 \\
7 & 7 & 8 \\
\end{array}
$
In contrast to step 1, we have a couple of exponents after the multiplication that occur twice. We are in modulo 2, meaning we remove those exponents that occur two times (2 is an even number). The exponents occurring twice are 5, 4, 3, 2, and 1.
Which leaves the exponents 8, 7, 6, and 0.
Create a polynomial with the exponents that remained after modulo 2:
$\mathtt{x^8 + x^7 + x^6 + x^0}$
That polynomial in binary and hex is
$\mathtt{111000001}_2$ and $\mathtt{1c1}_{16}$.
$\mathtt{1c1}_{16}$ is bigger than a byte, so we XOR $\mathtt{1c1}_{16}$ with $\mathtt{11b}_{16}$:
printf "%x\n" $((0x1c1 ^ 0x11b))
Which results in $\mathtt{da}_{16}$. We're done with step 2.
Steps 3 and 4
$\mathtt{5d} \cdot \mathtt{01}$ and $\mathtt{30} \cdot \mathtt{01}$ are easy because multiplication with 1 doesn't change the byte: We'll use $\mathtt{5d}$ and $\mathtt{30}$ for the final step where we XOR the four bytes from each step.
Step 5
In the final step, we add the bytes calculated in the four previous steps:
$$\mathtt{b3} + \mathtt{da} + \mathtt{5d} + \mathtt{30}$$
Adding bytes in $GF(2^8)$ means XORing them:
printf "%x\n" $((0xb3 ^ 0xda ^ 0x5d ^ 0x30))
We get $\mathtt{04}_{16}$ which is the result we expected:
$$(\mathtt{d4} \cdot \mathtt{02}) + (\mathtt{bf} \cdot \mathtt{03}) + (\mathtt{5d} \cdot \mathtt{01}) + (\mathtt{30} \cdot \mathtt{01}) = \mathtt{04}$$
Epilogue
If you're wondering why many implementations of MixColumns, such as the ones given on Wikipedia, use $\mathtt{1b}_{16}$ for reduction instead of $\mathtt{11b}_{16}$, this answer of mine might help you.
