clear
close all

m=1;
f=1;
xi=0.00;

w=2*pi*f;
k=m*w^2;
c=2*m*w*xi;
T=1/f;

P0 = 0;      % amplitude de l'excitation
f0 = 0.9;    % fréquence de l'excitation

x0=1;        % position initale
v0=0;        % vitesse initiale

Dt=0.01;     % pas de temps

N=round(15/Dt);
t=(0:N)*Dt;
p= P0 * sin(2*pi*f0*t);

%% Différence centrale 1

q(1)=x0;
v(1)=v0;
for i=1:N
    q(i+1)=q(i)+Dt*v(i)+Dt^2/2/m*(p(i)-c*v(i)-k*q(i));
    v(i+1)=2*(q(i+1)-q(i))/Dt-v(i);
end

showit = 1;  % utiliser 0 ou 1 pour le cacher ou le montrer
plot(t,showit*q); hold on

%% Différence centrale 1 - Version FF

q(1)=x0;
v(1)=v0;
for i=1:N
    q(i+1)=q(i)+Dt*v(i)+Dt^2/2/m*(p(i)-c*v(i)-k*q(i));
    v(i+1)=(q(i+1)-q(i))/Dt;
    %v(i+1)=v(i)+Dt/2/m*(p(i)-c*v(i)-k*x(i));
end

showit = 1;  % utiliser 0 ou 1 pour le cacher ou le montrer
plot(t,showit*q); hold on

%% Différence centrale 2
clear x

x_=x0+Dt*v0;
q(1)=x_;
q(2)=x0;

for i=2:N
    q(i+1)= (p(i)+(2*m/Dt^2-k)*q(i)-(m/Dt^2-c/2/Dt)*q(i-1)) / (m/Dt^2+c/2/Dt);
end

showit = 1;  % utiliser 0 ou 1 pour le cacher ou le montrer
plot(t,showit*q); hold on
legend('option 1','option 1-FF','option 2')

%ylim([-3 3])