In this guide to Macaulay2 I found an interesting way to compute intersection multiplicity. On page 61 the authors give an explicit case of the following argument:
We want to compute the intersection multiplicity at the origin of $n$ hypersurfaces in $\mathbb{A}^n_k$, with $k$ a field. Let $R=k[\mathbb{A}^n_k]$, $f_i$ be the equations of the hypersurfaces, $I=(f_1,\ldots,f_n)$, and let $m$ be the maximal ideal of the origin. Then this amounts to computing the dimension $\dim_k R_m/IR_m,$ which is the same as the dimension $\dim_k R/\varphi^{-1}(IR_m),$ where $\varphi$ is the natural localization map $\varphi: R\to R_m$.
Now the claim is that this ideal $\varphi^{-1}(IR_m)$ can be expressed as a quotient ideal $(I:(I:m^\infty)),$ where $(I:J)=\{x\in R\ |\ xJ\subset I\}$ and $(I:m^\infty)=\bigcup_{N=1}^{\infty}(I:m^N)=\{x\in R\ |\ \exists N:\ xm^N\subset I\}$. This is where I got stuck. It is a general fact that $\varphi^{-1}(IR_m)=\bigcup_{s\notin m}(I:s)$, but how do I equate it with the quotient ideal above? After this is done, one can compute the multiplicity as the degree of $(I:(I:m^\infty)).$
If so, then what additional assumption makes it true in the particular case covered in the book, where $n=3$, $I=(x^5 + y^3 + z^3,x^3 + y^5 + z^3,x^3 + y^3 + z^5)$?
I got it. The main idea is that $V(I:J)=\text{cl}(V(I)\setminus V(J))$: Closure of difference of varieties is ideal quotient.
Once we know that, we find that $V(I:m^\infty)$ consists of all components of $V(I)$ not containing the origin and $V(I:(I:m^\infty))$ is the origin with the correct scheme structure, and then the intersection multiplicity is just the degree of this 0-dimensional scheme.