I am answering my own question because I found the error. It was something in the code. I'm leaving the results here if anyone finds it useful or maybe wants to have a bit of fun with it:
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import numpy as np
w = np.linspace(-10, 10, 1000)
w1 = 2 * np.pi
A = np.array([-10, -5, -2, -1, -0.5])
B = -A
domain1 = np.linspace(A[0], B[0], 1000)
domain2 = np.linspace(A[1], B[1], 1000)
domain3 = np.linspace(A[2], B[2], 1000)
domain4 = np.linspace(A[3], B[3], 1000)
domain5 = np.linspace(A[4], B[4], 1000)
func1 = np.sin(w1 * domain1)
func2 = np.sin(w1 * domain2)
func3 = np.sin(w1 * domain3)
func4 = np.sin(w1 * domain4)
func5 = np.sin(w1 * domain5)
im1 = (np.sin((w1 - w) * A[0]) - np.sin((w1 - w) * B[0])) / (2 * (w1 - w)) + (
np.sin(-(w1 + w) * A[0]) - np.sin(-(w1 + w) * B[0])) / (2 * (w1 + w))
im2 = (np.sin((w1 - w) * A[1]) - np.sin((w1 - w) * B[1])) / (2 * (w1 - w)) + (
np.sin(-(w1 + w) * A[1]) - np.sin(-(w1 + w) * B[1])) / (2 * (w1 + w))
im3 = (np.sin((w1 - w) * A[2]) - np.sin((w1 - w) * B[2])) / (2 * (w1 - w)) + (
np.sin(-(w1 + w) * A[2]) - np.sin(-(w1 + w) * B[2])) / (2 * (w1 + w))
im4 = (np.sin((w1 - w) * A[3]) - np.sin((w1 - w) * B[3])) / (2 * (w1 - w)) + (
np.sin(-(w1 + w) * A[3]) - np.sin(-(w1 + w) * B[3])) / (2 * (w1 + w))
im5 = (np.sin((w1 - w) * A[4]) - np.sin((w1 - w) * B[4])) / (2 * (w1 - w)) + (
np.sin(-(w1 + w) * A[4]) - np.sin(-(w1 + w) * B[4])) / (2 * (w1 + w))
re1 = (np.cos((w1 - w) * A[0]) - np.cos((w1 - w) * B[0])) / (2 * (w1 - w)) + (
np.cos(-(w1 + w) * A[0]) - np.cos(-(w1 + w) * B[0])) / (2 * (w1 + w))
re2 = (np.cos((w1 - w) * A[1]) - np.cos((w1 - w) * B[1])) / (2 * (w1 - w)) + (
np.cos(-(w1 + w) * A[1]) - np.cos(-(w1 + w) * B[1])) / (2 * (w1 + w))
re3 = (np.cos((w1 - w) * A[2]) - np.cos((w1 - w) * B[2])) / (2 * (w1 - w)) + (
np.cos(-(w1 + w) * A[2]) - np.cos(-(w1 + w) * B[2])) / (2 * (w1 + w))
re4 = (np.cos((w1 - w) * A[3]) - np.cos((w1 - w) * B[3])) / (2 * (w1 - w)) + (
np.cos(-(w1 + w) * A[3]) - np.cos(-(w1 + w) * B[3])) / (2 * (w1 + w))
re5 = (np.cos((w1 - w) * A[4]) - np.cos((w1 - w) * B[4])) / (2 * (w1 - w)) + (
np.cos(-(w1 + w) * A[4]) - np.cos(-(w1 + w) * B[4])) / (2 * (w1 + w))
mod1 = np.sqrt(im1 ** 2 + re1 ** 2)
mod2 = np.sqrt(im2 ** 2 + re2 ** 2)
mod3 = np.sqrt(im3 ** 2 + re3 ** 2)
mod4 = np.sqrt(im4 ** 2 + re4 ** 2)
mod5 = np.sqrt(im5 ** 2 + re5 ** 2)
phase1 = np.arctan(im1 / re1)
phase2 = np.arctan(im2 / re2)
phase3 = np.arctan(im3 / re3)
phase4 = np.arctan(im4 / re4)
phase5 = np.arctan(im5 / re5)
fig = plt.figure()
grid = gridspec.GridSpec(nrows=3, ncols=2, figure=fig)
ax1 = fig.add_subplot(grid[0,:])
ax2 = fig.add_subplot(grid[1,0])
ax3 = fig.add_subplot(grid[1,1])
ax4 = fig.add_subplot(grid[2,0])
ax5 = fig.add_subplot(grid[2,1])
ax1.plot(w, mod1, label=r'$x\in [-10;10]$')
ax1.set_title(r'$|X(j\omega)|$')
ax1.grid(True)
ax1.vlines(w1, 0, 1, linestyles="dashed")
ax1.vlines(-w1, 0, 1, linestyles="dashed")
ax1.set_xlabel(r'$\omega$')
ax1.legend()
ax2.plot(w, mod2, label=r'$x\in [-5;5]$')
ax2.set_title(r'$|X(j\omega)|$')
ax2.grid(True)
ax2.vlines(w1, 0, 1, linestyles="dashed")
ax2.vlines(-w1, 0, 1, linestyles="dashed")
ax2.set_xlabel(r'$\omega$')
ax2.legend()
ax3.plot(w, mod3, label=r'$x\in [-2;2]$')
ax3.set_title(r'$|X(j\omega)|$')
ax3.grid(True)
ax3.vlines(w1, 0, 1, linestyles="dashed")
ax3.vlines(-w1, 0, 1, linestyles="dashed")
ax3.set_xlabel(r'$\omega$')
ax3.legend()
ax4.plot(w, mod4, label=r'$x\in [-1;1]$')
ax4.set_title(r'$|X(j\omega)|$')
ax4.grid(True)
ax4.vlines(w1, 0, 1, linestyles="dashed")
ax4.vlines(-w1, 0, 1, linestyles="dashed")
ax4.set_xlabel(r'$\omega$')
ax4.legend()
ax5.plot(w, mod5, label=r'$x\in [-0.5;0.5]$')
ax5.set_title(r'$|X(j\omega)|$')
ax5.grid(True)
ax5.vlines(w1, 0, 1, linestyles="dashed")
ax5.vlines(-w1, 0, 1, linestyles="dashed")
ax5.set_xlabel(r'$\omega$')
ax5.legend()
plt.tight_layout()
plt.show()
All this work follows from 3Blue1Brown video linking fourier transform and the quantum uncertainty principle. In the figure the legend tells us the time domain of x(t). Notice that the wider the time domain the more acute the maximums of the fourier transform get. That is because, if we are certain of where a particle is then its wave function will collapse to a dirac (I think so at least) and so its fourier transform will be a constant, that is, if we know exactly where a particle is, we have no idea what its frequency, and thus energy, is, it could be everywhere with the same likelihood in the frequency domain. As for the other limit, if we have no idea where a particle is, than its wave function spreads to infinity and the fourier transform converges to two diracs (positive and negative frequencies) giving us the exact frequency of the particle! Therefore, the the more certain we are of where a particle is, the thiner it's time domain has to be and its frequency domain is ever larger. We will never know for certain where and what its energy are at the same time. The fourier transform explains this concept beautifully.
In mathematics, there are many Fourier transforms. For every locally compact abelian group, there is a corresponding Fourier transform. Two of these Fourier transforms can be applied to a function $$f:[t_1,t_2]\rightarrow\mathbb{C}.$$
1) You can apply the Fourier transform $\mathscr{F}_{\mathbb{T}}$, where $\mathbb{T}:=\mathbb{R}/(t_2-t_1)\mathbb{Z}$ is a torus group. An application of $\mathscr{F}_{\mathbb{T}}$ corresponds to the action of extending $f$ periodically and then expanding it into a Fourier series. If you choose $\mathscr{F}_{\mathbb{T}}$ as your Fourier transform, $\mathscr{F}_{\mathbb{T}}(f)$ will be unique.
2) You can also choose to apply the Fourier transform $\mathscr{F}_{\mathbb{R}}$, but then you have to extend $f$ from $[t_1,t_2]$ to $\mathbb{R}$ beforehand. As you pointed out yourself, there is no unique way of doing this. Consequently, there is no unique way to define $\mathscr{F}_{\mathbb{R}}(f)$.
Addition (Discrete Fourier Transform):
If you consider the finite abelian group $\mathbb{Z}_N:=\mathbb{Z}/N\mathbb{Z}$, the corresponding Fourier transform $\mathscr{F}_{\mathbb{Z}_N}$ is called distrete Fourier Transform. If you only have $N$ sampling points $f(x_0),f(x_1),\ldots,f(x_{N-1})$, you may think of it as
a function on $\mathbb{Z}_N$ and consequently apply the Fourier transform
$\mathscr{F}_{\mathbb{Z}_N}$ to it.
I will not define the Fourier transforms $\mathscr{F}_{\mathbb{T}}$, $\mathscr{F}_{\mathbb{R}}$, $\mathscr{F}_{\mathbb{Z}_N}$ here (you can easily find good sources elsewhere), but simply emphasize that there is more than one Fourier transform, which may resolve some of the perplexity behind the question in the OP.
Best Answer
Let $f(x)$ be defined as
$$f(x)=\begin{cases}\cos(x)&|x|\le \pi/2\\\\0&,|x|\ge\pi/2\end{cases}$$
We can write $f(x)=\cos(x)\text{rect}(x/\pi)$ in terms of the rectangular function
$$\text{rect}(x)=\begin{cases} 1&, |x|<1/2\\\\0&,|x|>1/2 \end{cases}$$
Then, we have in distribution
$$f'(x)=-\sin(x)\text{rect}(x/\pi) $$
and
$$f''(x)=-\cos(x)\text{rect}(x/\pi)+\delta(x-\pi/2)+\delta(x+\pi/2)$$
So, the Fourier Transform of $f''(x)$ is
$$\begin{align} \mathscr{F}\{f''\}(\omega)&=\mathscr{F}\{-f\}+e^{-i\omega \pi/2}+e^{i\omega \pi/2}\\\\ &=-\frac{2\cos(\omega \pi/2)}{1-\omega^2}+2\cos(\omega \pi/2)\\\\ &=-\frac{2\omega^2 \cos(\omega \pi/2)}{1-\omega^2}\\\\ &=-\omega^2 \mathscr{F}\{f\}(\omega) \end{align}$$
as expected!