Replacing your line

[f,F2_x_f]=fourier(t_real,y,'sinus');

with

L=numel(t_real);
Fs=1/step_t;
f = Fs*(0:(L/2))/L;
F_x_f = fft(y,numel(f));
P2 = abs(F_x_f/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

I obtain a slightly different phase

zoomed in

I also moved the head line of graph 4 from Y axis to graph top,

..
% plot the phase spectrum
figure(3)
ax3=gca
plot(ax3,f,angle(F_x_f))
grid on
xlim([0,20])
xlabel('f[Hz]'); ylabel('[°]')
title('Phase spectrum of the force in °')
..

so all the readers of your question can comfortably read it.