%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% oscylator harmoniczny x'' + 2*p*x' + w0^2*x = F(t)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


p  = 0;		% tłumienie
w0 = 6;		% częstość własna
A  = 3;		% amplituda wymuszenia
w  = 4;		% częstość wymuszenia
T  = 100;	% czas obserwacji

F = @(t) A*sin(w*t); % zewnętrzne wymuszenie

% zamieniamy na układ równań:
%
% x' = v
% v' = F(t) - 2*p*v - w0^2*x
%
% X = [X(1), X(2)] = [x,v]
f = @(X,t) [X(2); F(t) - 2*p*X(2) - (w0^2)*X(1)];

t = linspace(0,T,3000);
X = lsode(f,[0;0],t);
% wykres rozwiązania
plot(t,X); grid on; legend('x(t)', 'v(t)'); pause(1.5);
% wykres fazowy
plot(X(:,1),X(:,2)); xlabel('x(t)'); ylabel('y(t)'); pause(1.5);


