Funcion para el calculo de la respuesta dinamica utilizando el metodo de interpolacion de la exitacion (Time Stepping Method).
[xt,vt,coef,Tabla] = Time_Stepping(m,k,xi,t,P,delta_t,xo,vo)
A partir de la solucion de ecuaciones de movimiento en el intervalo de tiempo, tenemos:
u(i+1)=Aui+Bvi+CPi+DP(i+1)
v(i+1)=A1ui+B1vi+C1Pi+D1P(i+1)
[Tabla] = [ti pi ui vi u(i+1) v(i+1)]
function [xt,vt,coef,Tabla] = Time_Stepping(m,k,xi,t,P,delta_t,xo,vo)

wn=(k/m)^0.5;
wd=wn*(1-xi^2)^0.5;

% Calculo de los coeficientes:
% Si el delta_t es constante los coeficientes deben ser calculados una sola
% vez.

A=exp(-xi*wn*delta_t)*((xi/(1-xi^2)^0.5)*(sin(wd*delta_t))+(cos(wd*delta_t)));
B=exp(-xi*wn*delta_t)*((1/wd)*sin(wd*delta_t));
C=(1/k)*((2*xi/(wn*delta_t))+exp(-xi*wn*delta_t)*((((1-2*xi^2)/(wd*delta_t))-(xi/(1-xi^2)^0.5))*(sin(wd*delta_t))-(1+2*xi/(wn*delta_t))*(cos(wd*delta_t))));
D=(1/k)*(1-(2*xi/(wn*delta_t))+exp(-xi*wn*delta_t)*(((2*xi^2-1)/(wd*delta_t))*(sin(wd*delta_t))+(2*xi/(wn*delta_t))*(cos(wd*delta_t))));

A1=-exp(-xi*wn*delta_t)*((wn/(1-xi^2)^0.5)*(sin(wd*delta_t)));
B1=exp(-xi*wn*delta_t)*(cos(wd*delta_t)-(xi/(1-xi^2)^0.5)*(sin(wd*delta_t)));
C1=(1/k)*((-1/delta_t)+exp(-xi*wn*delta_t)*(((wn/(1-xi^2)^0.5)+(xi/(delta_t*(1-xi^2)^0.5)))*(sin(wd*delta_t))+(1/delta_t)*(cos(wd*delta_t))));
D1=(1/(k*delta_t))*(1-exp(-xi*wn*delta_t)*((xi/(1-xi^2)^0.5)*(sin(wd*delta_t))+(cos(wd*delta_t))));

coef=[A B C D
      A1 B1 C1 D1]';

% Definimos los vectores de posicion y velocidad
%             1  2  3  4     5      6    
% [Tabla] = [ti pi ui vi u(i+1) v(i+1)]
% u(i+1)=Aui+Bvi+CPi+DP(i+1)
% v(i+1)=A1ui+B1vi+C1Pi+D1P(i+1)

Tabla=zeros(length(t),8);

Tabla(:,1)=t';
Tabla(:,2)=P';
Tabla(1,3)=xo;
Tabla(1,4)=vo;

for i=1:length(t)-1

    Tabla(i,5)=A*Tabla(i,3)+B*Tabla(i,4)+C*Tabla(i,2)+D*Tabla(i+1,2);
    Tabla(i+1,3)=Tabla(i,5);

    Tabla(i,6)=A1*Tabla(i,3)+B1*Tabla(i,4)+C1*Tabla(i,2)+D1*Tabla(i+1,2);
    Tabla(i+1,4)=Tabla(i,6);

end

xt=Tabla(:,3);
vt=Tabla(:,4);

end