Plot is not showing anything.

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

for t = 0:h:1

xt=X0*exp((-S)*wn*t)*cos((wd*t)-Phi0)+X*cos((w*t)-Phi);

t = t+h; %time

end

figure(1)

plot(xt,t);

xlabel('time');

ylabel('Complete Solution');

title('Second Order ODE');

hold on

Answer 1

[Les Beckham]

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)