Replace your loop with
h= 0.1; %step size for loop
t = 0:h:1;
for i = 1:numel(t)
xt(i) = X0*exp((-S)*wn*t(i))*cos((wd*t(i))-Phi0) + X*cos((w*t(i))-Phi);
end
so that you save xt for every time sample.
k= 4000; %N/m
m= 10; %kg
w= 10; %forcing frequency (rad/sec)
F0= 100; %applied force (N)
c= 20; %(N*s)/m Damping constant
x0=0.01; %meters
x0dot=0; %meters/sec
wn= sqrt(k/m); %natural frequency (rad/sec)
cc= 2*(sqrt(k*m)); %critical damping constant
S= c/cc; %Damping ratio
wd= sqrt(1-((S^2)*wn)); %damping frequency
r= w/wn; %frequency ratio
Sdef= F0/k; %Static Deformation
X=(Sdef)/(sqrt(((1-(r^2))^2)+((2*S*r)^2))); %Amplitude
Phi= atan((2*S*r)/(1-(r^2))); %Phase Angle
X0=sqrt(((x0-X*cos(Phi))^2)+((1/(wd^2))*(((S*wn*x0)+(x0dot)-(S*wn*X*cos(Phi))-(w*X*sin(Phi)))^2)));
Phi0= atan(((S*wn*x0)+(x0dot)-(S*wn*X*cos(Phi))-(w*X*sin(Phi)))/((wd)*(x0-X*cos(Phi))));
h= 0.1; %step size for loop
t = 0:h:1;
for i = 1:numel(t)
xt(i) = X0*exp((-S)*wn*t(i))*cos((wd*t(i))-Phi0) + X*cos((w*t(i))-Phi);
end
figure(1)
plot(xt,t);
xlabel('time');
ylabel('Complete Solution');
title('Second Order ODE');
Note that you don't even need the loop. Matlab will do this as a "vectorized" operation all at once.
xt = X0 * exp((-S)*wn*t) .* cos((wd*t)-Phi0) + X*cos((w*t)-Phi);
plot(xt, t)
