clc 
clear all 

%% params 
m = 1;
k = 4;
c = 0.1;

fsamp = 512;
Tmax = 0.01;
fpuls = 20;

T = linspace(0, Tmax, Tmax*fsamp+1)';
fex = sin(2*pi*fpuls*T);

%% LSIM
A = [0 1;-k/m -c/m];
B = [0; 1/m];
C = eye(2);

sys = ss(A,B,C,0);
[Y, T] = lsim(sys, fex, T, [0;0]);
[To, Yo] = ode45(@(t,y) A*y+B*sin(2*pi*fpuls*t), T, [0;0]);


figure(5)
clf()
plot(T, Y(:,2), '.-', To, Yo(:,2), 'o-'); grid on
legend('lsim', 'ode45', 'Location', 'northwest')
xlabel('Time (s)')
ylabel('Response (m/s)')
title('Octave')
% title('Matlab')
