Funcion para el calculo de la respuesta dinamica utilizando el metodo de diferencias finitas (Time Stepping Method).
[xt,vt,at,Tabla,Coeficientes] = Diferencia_Central(m,k,xi,t,P,delta_t,x0,v0)
function [xt,vt,at,Tabla,Coeficientes,Check] = Diferencia_Central(m,k,xi,t,P,delta_t,x0,v0)

% Calculos iniciales

wn=(k/m)^0.5;

c=xi*(2*m*wn);
P0=P(1,1);
a0=(P0-c*v0-k*x0)/(m);
x_menos1=x0-v0*delta_t+((delta_t)^2/2)*a0;
k_hat=(m/delta_t^2)+(c/(2*delta_t));
a=(m/delta_t^2)-(c/(2*delta_t));
b=k-(2*m/delta_t^2);

Coeficientes=[c P0 a0 x_menos1 k_hat a b];

% Debemos resolver la ecuacion
% k_hat*x(i+1)=Pi-a*x(i-1)-b*x(i)

P=[P0 P]';
t=[-delta_t t]';

% Tabla = [t P ut vi ai]
Tabla=zeros(length(t),2);
Tabla(1,3)=x_menos1;
Tabla(2,3)=x0;
Tabla(:,1)= P;
Tabla(:,2)= t;

for i=3:length(t)
    % En este caso Tabla(i,3) = x(i+1)
    % Por lo tanto Tabla(i-1,3) = x(i)
    %              Tabla(i-2,3) = x(i-1)
    Tabla(i,3)=(Tabla(i-1,1)-a*Tabla(i-2,3)-b*Tabla(i-1,3))/(k_hat);
    Tabla(i,4)=(Tabla(i,3)-Tabla(i-2,3))/(2*delta_t);
    Tabla(i,4)=0;
end

% Eliminamos los valores del tiempo menos 1

xt=Tabla(2:end,3);
vt=Tabla(2:end,4);
at=Tabla(2:end,4);

Tn=2*pi/wn;

if delta_t/Tn<1/pi
    Check="OK";
else
    Check="No cumple";
end

end