Say I have the following operator:
$$\hat { L } =\hbar { \sum_{ \sigma ,l,p } { l } \int_{ 0 }^{ \infty }\!{ \mathrm{d}{ k }_{ 0 }\,\hat { { { a }}}_{ \sigma ,l,p }^{ \dagger } } } \left({ k }_{ 0 }\right)\hat { { a }}_{ \sigma ,l,p } \left(k _{ 0 }\right)$$
How would I write this in matrix form?
The question is a little vague but you first need an orthonormal basis, i.e. you need a set of vectors $\vert \phi_k\rangle$ for which you can compute $\langle \phi_j\vert \hat L\vert \phi_k\rangle = c_{jk}$. Then it’s a matter of writing $$ \hat L=\sum_{j’,k’} c_{j’k’} \vert \phi_{j’}\rangle\langle \phi_{k’}\vert \tag{1} $$ so that $$ \langle \phi_j\vert \hat L\vert\phi_k\rangle = \sum_{j’,k’} c_{j’k’}\langle \phi_j\vert\phi_{j’}\rangle \langle \phi_{k’}\vert \phi_k\rangle= c_{jk} $$ follows from (1) by orthonormality $\langle \phi_{j’}\vert\phi_j\rangle=\delta_{j’j}$. The $c_{jk}$ can then immediately be inserted into a matrix at position $(j,k)$.
Altough this is more mathematics than physics, I'll leave this comment:
You need to know a basis for the space. For eaxmple, in a 2-D space, you might choose the $\left( \begin{array}{} 1 \\ 0 \end{array} \right) ; \left( \begin{array}{} 0 \\ 1 \end{array} \right) $ basis.
Finally, you'd have to know how the operator transforms the initial-basis vectors. That means, the results of $L\left( \begin{array}{} 1 \\ 0 \end{array} \right) ;\quad L\left( \begin{array}{} 0 \\ 1 \end{array} \right) $
Then, it's basic algebra that the columns of the matrix are the the components of those results.
The first question is: What's the vector space? You need to know where you are in order to find a basis. Then you should check that you're working with functions, which are infinite dimensional vector spaces. You cannot build an infinite matrix!
So trying to find such matrix is absurd You just work with the differential operator.
An operator is a linear function $f \in H \to Y$ where $H$ and $Y$ are vector spaces. For such a function there are representations that resemble matrices. For example, it is common in physics that $H=Y$ is a Hilbert space over the complex numbers and that is has a Schauder basis, in the sense that is has a countable set $(y_n)_{n=1}^\infty$ such that for any element $x \in H$ there is a sequence of of complex numbers $(a_n)_{n=1}^\infty$ such that $x=\sum_{n=1}^\infty a_ny_n$, the sum converging in the Hilber space norm. In this case, $f(x)=\sum_{n=1}^\infty a_nf(y_n)$ and also there is a sequence $(b_{n,m})_{m=1}^\infty$ such that $f(y_n)=\sum_{m=1}^\infty b_{n,m}y_m$ so
$f(x)=\sum_{n=1}^\infty\sum_{m=1}^\infty a_nb_{n,m}y_m$
which looks like how a matrix operates. You could represent that as a countably infinite matrix $(b_{n,m})_{n=1,m=1}^\infty$ being multiplied by other such matrices in the usual way and you can check that composition of operators is respected by this representation.
However this is but an example. Another one would consider bigger sets of "orthonormal" vectors and you would have integrals instead of summations (and extra trouble?). You might be using another "sense" of a basis and you would have to go through a similar procedure to build the representation. In general, once you have the representation and have chosen your "basis" the "matrix elements" would be given by the projections of $f(y_i)$ on $y_j$, $P_{y_j}f(y_i)$ where y stands for a "basis" element and its index is your way of identifying each of them (a counting set, for example the real numbers).
