Analisis Ecuacion de movimiento

Introduccion
La tarea consiste en resolver la ecuación de movimiento sujeta a cargas periódicas que pueden ser escritas en forma de sumas de senos y cosenos mediante la serie de Fourier.
MarcoTeórico
La ecuación a resolver es la siguiente:
mv ̈+cv ̇+kv=p(t),con p(t)=a_0+∑_(n=1)^∞▒〖a_n cos⁡(2πn/T_p t) 〗+b_n sin⁡(2πn/T_p t)
La función p(t) es la serie de Fourier de un pulso P(t), queviene dado por el siguiente gráfico:

Gráfico 1: Carga P(t).
La solución a la ecuación anterior viene dada por la siguiente expresión:
v(t)=1/k [a_o+∑_(n=1)^∞▒〖[1/((1-(ω ̅_n/ω)^2 )^2+2β(ω ̅_n/ω)^2)]×〗 {[2β ω ̅_n/ω a_n+b_n (1-(ω ̅_n/ω)^2 )]sen(ω ̅_n t)+[-2β ω ̅_n/ω b_n+a_n (1-(ω ̅_n/ω)^2 )]cos(ω ̅_n t)}]
Ecuación 1: Posición en función del tiempo.
En donde:
a_o=1/T_p ∫_0^Tp▒P(t)dt; a_n=2/T_p∫_0^Tp▒〖P(t)cos⁡(2πt/T_p )dt〗; b_n=2/T_p ∫_0^Tp▒〖P(t) sin(2πt/T_p )dt; 〗
〖y ω ̅〗_n=2πn/T_p
Ecuaciones 2, 3, 4, 5, Parámetros de la ecuación 1.
Derivando la ecuación 1 se obtiene la velocidad:
v̇(t)=1/k [∑_(n=1)^∞▒〖[1/((1-(ω ̅_n/ω)^2 )^2+2β(ω ̅_n/ω)^2 )]×〗 {[2β ω ̅_n/ω a_n+b_n (1-(ω ̅_n/ω)^2 )] ω ̅_n cos(ω ̅_n t)-[-2β ω ̅_n/ω b_n+a_n (1-(ω ̅_n/ω)^2 )] ω ̅_n sin(ω ̅_n t)}]
Ecuación 6:Posición en función del tiempo.
v ̈(t)=1/k [∑_(n=1)^∞▒〖[1/((1-(ω ̅_n/ω)^2 )^2+2β(ω ̅_n/ω)^2 )]×〗 {-[2β ω ̅_n/ω a_n+b_n (1-(ω ̅_n/ω)^2 )] ω ̅_n^2 sin(ω ̅_n t)-[-2β ω ̅_n/ω b_n+a_n (1-(ω ̅_n/ω)^2 )] ω ̅_n^2cos(ω ̅_n t)}]
Ecuación 7: Aceleración en función del tiempo.
La resolución de la ecuación 1 se realizó en MATLAB, se realizó un ciclo for, para ir sumando los términos de la suma que va desde cerohasta infinito (este infinito se considerara como 20, que es una buena aproximación, al mismo tiempo este ciclo for va calculando los coeficientes a_n y b_n de la serie de Fourier utilizando lafunción trapz, que calcula integrales numéricamente.
Las consideraciones para resolver el problema son las siguientes:
Los datos conocidos son los siguientes:
m=1
β=3%
T=0.4 [s]
Las...