In practice, you don't. You could do it. An output of a layer of a neural net is just a bunch of linear combinations of the input followed by a (usually non-linear) function application (a sigmoid or, nowadays ReLU). Then you just do this again for each layer. You can easily write out what this equation must be. You could symbolically differentiate that, but the equation is massive. You can simplify somewhat and recognize that the output is a composition of functions that I described above and so you can write its derivative as multiple applications of the chain rule. You can then see what you'd need to do to calculate the factors of the resulting expression layer by layer. This leads to the backpropagation algorithm. Then you talk to the right kind of mathematician (like a numerical analyst) and they tell you that backpropagation is just a special case of reverse mode automatic differentiation. Automatic differentiation is a readily implementable technique that allows you to turn a fairly arbitrary program that calculates a mathematical function, into a program that calculates that function and its derivative. In other words, the entire backpropagation idea of neural nets can be reduced to: 1) write an program that calculates the value of the neural net, 2) apply automatic differentiation to it to get its derivative, 3) do the obvious gradient descent thing (i.e. move a bit in the direction the gradient suggests), 4) repeat.
Automatic differentiation doesn't require or produce an explicit expression for the derivative, nor does it approximate the derivative. The basic idea (for forward mode AD at least) is straightforward. Write $Df$ for the derivative of $f$ with respect to its argument. Let's say we wanted to know what the derivative of $f+g$ is at $x$, i.e. $D(f+g)(x)$. We know from basic calculus that this is $(Df+Dg)(x)=Df(x)+Dg(x)$. So if the program implementing $f$ when evaluated at $x$ produced not only the value $f(x)$ but also the (single numerical value!) $Df(x)$ and similarly for $g$, we could calculate $D(f+g)(x)$ by simply adding those two extra outputs. Similarly, for $D(fg)(x)=f(x)Dg(x)+g(x)Df(x)$ we use both the "normal" outputs of $f$ and $g$ and the "extra" derivative outputs and easily calculate the the "extra" derivative output of the product of the functions. We can do similar things for various other combinations of functions (scaling, composition, exponentiation) and easily implement primitive operations (like $\sin$ and $\cos$) to produce these "extra" derivative values. With operator overloading, type classes, or program rewriting, you can just work in terms of the "normal" values and automatically, in parallel, the derivatives will also be calculated. Reverse mode AD is a little more complicated but the end experience is much the same.
So to reiterate, backpropagation is an algorithm that can be automatically derived and generated. You can implement forward mode automatic differentiation in Haskell, for example, in a few dozen lines of code, most of which are just writing out the derivatives of primitive operations. However, as I mentioned, backpropagation is reverse mode automatic differentiation which is harder to implement. It can still be done as a library in Haskell, but most implementations of reverse mode AD work as program transformations.