Funcion para el calculo de la respuesta dinamica utilizando el metodo de interpolacion de la exitacion (Time Stepping Method).
[xt,vt,at,coef,Tabla] = Newmark(m,k,xi,t,P,delta_t,x0,v0,beta,gamma)
Aceleracion constante - beta=1/4 y gamma=1/2
Aceleracion lineal - beta=1/6 y gamma=1/2
function [xt,vt,at,coef,Tabla,Check] = Newmark(m,k,xi,t,P,delta_t,x0,v0,beta,gamma)

wn=(k/m)^0.5;
Tn=2*pi/wn;

% Calculos iniciales

c=xi*(2*m*wn);
P0=P(1,1);
a0=(P0-c*v0-k*x0)/(m);
a1=((1)/(beta*delta_t^2))*m+((gamma)/(beta*delta_t))*c;
a2=(1/(beta*delta_t))*m+(gamma/beta-1)*c;
a3=(1/(2*beta)-1)*m+delta_t*(gamma/(2*beta)-1)*c;
k_hat=k+a1;

coef=[a1 a2 a3 k_hat]';

%           1  2  3  4  5    6
% Tabla = [ti Pi ui vi ai pi+1]
Tabla=zeros(length(t),6);
Tabla(:,1)=t;
Tabla(:,2)=P;
Tabla(1,3)=x0;
Tabla(1,4)=v0;
Tabla(1,5)=a0;

% Tenemos que resolver la ecuacion

for i=2:length(t)

p_hat_1=Tabla(i,2)+a1*Tabla(i-1,3)+a2*Tabla(i-1,4)+a3*Tabla(i-1,5);
Tabla(i,6)=p_hat_1;
Tabla(i,3)=p_hat_1/k_hat;
Tabla(i,4)=(gamma/(beta*delta_t))*(Tabla(i,3)-Tabla(i-1,3))+(1-gamma/beta)*Tabla(i-1,4)+delta_t*(1-gamma/(2*beta))*Tabla(i-1,5);
Tabla(i,5)=(1/(beta*delta_t^2))*(Tabla(i,3)-Tabla(i-1,3))-(1/(beta*delta_t))*Tabla(i-1,4)-(1/(2*beta)-1)*Tabla(i-1,5);
end

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

% Verificacion de estabilidad

    Lim=(1/(pi*sqrt(2)))*(1/(beta*sqrt(gamma-2*beta)));
    if delta_t/Tn < Lim
        Check="Ok";
    else
        Check="Not OK";
    end