fun = @(t, X) odefun(t, X, N, marray(m), carray(c, N), make_diagonal(X, k, gamma, N));
[t, X] = ode45(fun, tspan, X0);
plot(t/pi, X), grid on, xlabel('t/\pi'), ylabel('\bf{x}\rm(t)'), title('System responses')
function dX = odefun(t, X, N, M, C, K)
f = famp*sin(alpha*t)*cos(beta*t);
F = [f; zeros(175,1); f; zeros(100,1); f; zeros(74,1); f; zeros(175,1); f; zeros(71,1)];
ddx = M\(F - K*x - C*dx);
function K = karray(k, gamma, N)
K = @(x) make_diagonal(x, k, gamma, N);
function out = make_diagonal(x, k, gamma, N)
cg = circshift(gamma, -1);
d1 = -3.*ck.*cg.*cx.^2 - ck;
d2 = (k.*gamma + ck.*cg).*x.^2 + k + ck;
out = full(spdiags([d1 d2 d3], -1:1, N, N));