"Low pass" is a qualifier. I can recall maybe more then 5 low pass qualities but each bringing along a different phase response making it attractive to a particular application. Others may be able to recall even more.
What I want is not a ready-to-use algorithm that I can copypast. I want to understand what I'm doing and how it works.
This would make the question very broad. It is good however that you say that you have an understanding of convolution. So, if I had to put it in one phrase, all you do is pitching the right phases against each other. For more information, you might want to refer to this question and its answers.
Now, once you understand that interplay of bounces and phases, you have Finite Impulse Response (FIR) filters covered. This extends to Infinite Impulse Response (IIR) filters by noticing that the two branches of feed forward and feed back are essentially the sum of two convolutions back to back. This is more "visible" if you look at the direct form implementations.
But the idea is the same: pitch phases against each other.
How you do that? That would take a bit more space to explain but there are monumental works in this field that contain much more information.
However! All the documentation that I've found until now, including the Wikipedia pages, are directed to people that already have a solid background in signal processing (a lot of unexplained notations, function names, etc.). This is not my case
The notation is simply discrete sequences with the addition of the convolution operator ($*$). Within the context of DSP there are some standard signals used for specific purposes that sometimes are not explained. For example, the unit impulse or the step or the ramp and others. And also the standard notation used in transforms.
I would expect that you would pick these things up very quickly, especially if you follow one of those textbooks that are highly regarded.
I need someone to explain to me how a low pass filter works, how the Butterworth filter (or something else I don't care) works, how I can extract a nice algorithm from it with the correct coefficients. And all that like I'm a 5 years old.
Indeed, there is a little bit of contradiction here. You don't want an algorithm, you don't want to copy-paste but you want to extract an algorithm that gives you the coefficients.
Again, the "why" is answered by looking at the purpose for which the filter was designed for. For example, the Gaussian filter has "... no overshoot to a step function input while minimizing the rise and fall time...". So, it is the fastest non-ringing filter. That is its unique characteristic (notice here, still a low pass).
I appreciate that you need a starting point, so, one of the easiest methods to design an FIR filter is to design its frequency response in the frequency domain and obtain your impulse response, usually referred to as $h[n]$ (it is in the time domain), via the inverse Discrete Fourier Transform with the potential application of a window function to reduce the length of the filter gracefully (also referred to as the order of the filter).
The reason for wanting to reduce the number of coefficients is to reduce computational complexity. But you don't want to simply chuck away some coefficients, you want to do it in a controlled way so that the filter doesn't deviate too much from its specification.
One of the simplest expressions of this method is the sinc filter.
Which, in a few phrases, is explained as follows: Draw a square pulse in the frequency domain that starts from DC (0 Hz) and stops at the cut-off frequency of your low pass filter. Do an inverse DFT of that and this results to a sinc pulse. Use that sinc to filter your data. The longer the sinc, the sharper the cut off.
For more information about this, please see here, here, or here.
Here is a back of the envelope quick demo of this (in Octave but easily translatable to other platforms):
Fs=44100; % Sampling frequency in Hertz
N=44100; % Length of the impulse response (please see below) in samples
f=2000; % Cut off frequency of the filter
t=(-(N./2):(N./2))./Fs; % Time vector such that the sinc is always centered in the impulse response buffer
p=2*pi*t; % Phase vector (to save us multiplying by 2*pi everywhere)
h = sinc((fc./pi).*p); % Impulse response (we don't do IDFT because we know analytically its result and use the sinc function directly).
h = h./sum(h); % We do this to normalise the amplitude to 1. In other words, the filter will not amplify or attenuate the incoming signal by a constant factor.
semilogy(abs(fft(h,Fs)));grid on;xlabel("Frequency (Hz)"); ylabel("Amplitude");
The above, sets up a few variables and plots the frequency response of the filter which looks like this:

Which, as you can see, has a flat area in the beginning, that is your passband, followed by a very sharp roll-off region (the part that is almost vertical), followed by a stopband which shows a very high attenuation of frequencies beyond the 2kHz specification.
If you wanted to try this, pick a mono wav file (for convenience) that has the same sampling frequency, load it as a "signal" and then listen to the result of convolution. Like this:
% Given the h from above
s = wavread("myfile.wav");
o = conv(s,h);
plr = audioplayer(o,Fs);
play(plr)
(More information on the specifics of the above functions are available here).
Now, depending on the sound file, you might find that a cut off frequency of 2kHz isn't too audible so try f=200
above.
You might want to try different f
and N
just to see the effect they have on the filter.
This brushes over some key details (such as normalised frequency, transfer function and poles / zeros amongst others) but hopefully it will be enough to get you going.
Hope this helps