plot(t,c)
You are not storing into c in the workspace that you are doing the ode45 call, so it will have whatever value it had before the call.
function f = velocity(t,v)
f = 8*v/t
M = [4 5;5 6]
K = [23 45;67 54]
[X e] = polyeig(K,M)
F = [32+f;32+f]
wn = 3221.3;
zeta = 0.061;
wd = 3254.4;
Fn = X.'*F
Fnn = Fn(2,1)
g = -wd.^2*sin(wd*t) - 2*wd*zeta*wn*cos(wd*t) +zeta.^2*wn.^2*sin(wd*t)
a = Fnn.*(cos(-zeta*wd*t)+sin(-zeta*wd*t))/(8*wd.^2 +1i*2.12*wd+253)
c = a.*g
end
The only line there that affects the ode45 integration is the very first one, f = 8*v/t . All of the rest are wasted computation as far as ode45 is concerned.
I would suggest,
[t,v] = ode45(@velocity, [10 1], [22.5]);
c = arrayfun(@calc_c, t, v);
plot(t, c);
function f = velocity(t,v)
f = 8*v/t;
end
function c = calc_c(t,v)
f = 8*v/t;
M = [4 5;5 6];
K = [23 45;67 54];
[X e] = polyeig(K,M);
F = [32+f;32+f];
wn = 3221.3;
zeta = 0.061;
wd = 3254.4;
Fn = X.'*F;
Fnn = Fn(2,1);
g = -wd.^2*sin(wd*t) - 2*wd*zeta*wn*cos(wd*t) +zeta.^2*wn.^2*sin(wd*t);
a = Fnn.*(cos(-zeta*wd*t)+sin(-zeta*wd*t))/(8*wd.^2 +1i*2.12*wd+253);
c = a.*g;
end
