function [zs , alphas] = soluedo( RS,NPTs,dt,phi0s )
zs(1)=phi0s(1);
zs(2)=zs(1)+dt*phi0s(2);
alphas(1)=phi0s(3);
alphas(2)=alphas(1)+dt*phi0s(4);
X1=[zs(1),alphas(1)]';
X2=[zs(2),alphas(2)]';
%=============
% resolution du systeme edo
%===========
for i=3:NPTs
X3=2*X2-X1-dt^2*RS*X2;
zs(i)=X3(1);
alphas(i)=X3(2);
X1=X2;
X2=X3;
end