Hello
Following link helped me to solve my issue: https://www.12000.org/my_notes/on_scaling_factor_for_ftt_in_matlab/
I have divided the fft by the sampling time. Then I have also deleted the factor 2 at the magnitude and at the phase calculation. It is still not reasonable to me why, because I am plotting just one half (as I understood: frequencies from 0 to fs/2 and not -fs/2 to fs/2) of the fft. In order to keep the "signal energy" constant I should multiply by 2. Anyways, it works now without the "2".
clear all, clc, close all
%.........
% System
%.........
s=tf('s');
stoerung_G = 75 / (s + 25);
figure; bode(stoerung_G); grid;
stoerung_dt = 0.00001;
stoerung_t = [0 : stoerung_dt : 0.4]';
[stoerung_impulse] = impulse(stoerung_G, stoerung_t);
%.........
% Fourieranalyse
%.........
stoerung_y = stoerung_impulse;
stoerung_Fs = 1/stoerung_dt; % Hz
stoerung_L = size(stoerung_y,1);
stoerung_Y_ft = fft(stoerung_y)*stoerung_dt;
P2 = abs(stoerung_Y_ft);
P1 = P2(1:stoerung_L/2+1);
P1_log = 20 * log10(P1);
Phase2 = angle(stoerung_Y_ft);
Phase1 = Phase2(1:stoerung_L/2+1);
Phase1 = Phase1*180/pi();
stoerung_f = stoerung_Fs*(0:(stoerung_L/2))/stoerung_L;
stoerung_w = 2*pi()*stoerung_f;
figure
subplot(211)
semilogx(stoerung_w,P1_log), grid, xlim([1 3000]);
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('w (rad/s)')
ylabel('|Mag| (dB)')
subplot(212)
semilogx(stoerung_w,Phase1), grid, xlim([1 3000]);
xlabel('w (rad/s)')
ylabel('Phase (degree)')