M=1;   % Masa
k=1;   % Constante elastica del resorte
b=10;  % Coeficiente de roce
x1(1)=0; % posicion inicial de la masa
x2(1)=0; % velocidad inicial de la masa
T=0.01;  % paso de integracion numerica
Tfinal=30; % Tiempo final de integracion
t=[0:T:Tfinal]';
F=ones(size(t));
for n=2:length(t)
    x1(n)=T*x2(n-1)+x1(n-1);
    x2(n)=(T/M)*(F(n-1)-k*x1(n-1)-b*x2(n-1))+x2(n-1);
end
plot(t,x1,t,x2)
grid