Прохождение единичного скачка через колебательный контур

u(1) = 1; % единичный скачок

clc; clear all; close all;
R = 1;  % сопротивление, Ом
L = 1;  % индуктивность, Генри
C = 1e-3; % ёмкость, Фарады
Fs = 1000;  % частота дискретизации, Гц
T = 1;      % длительность сигнала, сек
t = transpose(0:1/Fs:20); % временная ось
Count = length(t);
Dirak = 0; %переключатель ф-ций
if Dirak
    % Дельта-ф-ция
    u = zeros(size(t));
    u(1000) = 1; 
else % Переходная ф-ция
    u = ones(size(t));
    u(1:500) = 0; 
    u(8000:20000) = 0; 
end
U = fft(u);
f = transpose(0:Count-1)*Fs/Count;
K = 1 ./ (1 - (2*pi*f).^2 * L * C + 1i * 2*pi*f * R * C); %частотная характеристика контура
Y = K .* U;%отклик в частотной области
y = ifft(Y); % напряжение на конденсаторе в временной области
V = real(y); % напряжение на конденсаторе
I = [0;diff(V)]*Fs; % ток в контуре
 
figure;
subplot(4, 2, 1);
plot(t, u);
title('Временная область')
ylabel('U(t)');
xlabel('t, s');
grid on;
 
subplot(4, 2, 2);
%heaviside1 = double(t >= 0);
plot(f, abs(U));
title('Частотная область')
ylabel('H_U(f)');
xlabel('f, Hz');
xlim([0 100]);
%ylim([0 1]);
grid on;
 
subplot(4, 2, 4);
%heaviside1 = double(t >= 0);
plot(f, abs(K));
%title('')
xlabel('f, Hz');
ylabel('К(f)');
xlim([0 100]);
grid on;
 
subplot(4, 2, 5);
plot(t, V);
ylabel('U_C(t)');
xlabel('t, s');
grid on;
 
subplot(4, 2, 6);
%heaviside1 = double(t >= 0);
plot(f, abs(Y));
xlabel('f, Hz');
ylabel('Y(f)');
xlim([0 100]);
grid on;
 
subplot(4, 2, 7);
plot(t, I);
ylabel('I_C(t)');
xlabel('t, s');
grid on;