I think the following produces somewhat more sensible results. I was confused for some time because you use x for both distance along the road and for vertical displacement. I've changed the latter to y. See if this does what you want:
k_s = 26400; %spring stiffness
m = 483; %Mass
f_n = sqrt(k_s/m)/(2*pi); %Natural frequency in Hz
%% Road profile
% spatial frequency (n0) cycles per meter
Omega0 = 0.1; %%%%conventional value of spatial frequency(n0)?
% psd ISO (used for formula 8)
Gd_0 = 32 * (10^-6);
% waveviness
w = 2;
% road length
L = 250;
%delta n
N = 100;
Omega_L = 0.004;
Omega_U = 4;
delta_n = 1/L; % delta_n = (Omega_U - Omega_L)/(N-1);
% spatial frequency band
Omega = Omega_L:delta_n:Omega_U;
%PSD of road
Gd = Gd_0.*(Omega./Omega0).^(-w);
% calculate amplitude using formula(8) in the article
%Amp = sqrt(2*Gd*delta_n); %%%from Eq. 7?
%calculate amplitude using simplified formula(9) in the article
k = 3; %%%upper limit A and lower limit B k=3?
%Amp = sqrt(delta_n) * (2^k) * (10^-3) * (Omega0./Omega);
Amp = sqrt(delta_n) * (2^k) * (10^-3) * (Omega0./Omega);
%random phases
Psi = 2*pi*rand(size(Omega));
% x abicsa from 0 to L
x1 = 0:250/(N-1):250;
h= zeros(size(x1));
%artificial random road profile
for iv=1:length(x1)
h(iv) = sum( Amp.*cos(2*pi*Omega*x1(iv) + Psi) );
end
hx = [x1' h'];
%% ode45
y0 = [0,0];
[t, y] = ode45(@f,x1,y0,[],hx);
%% plot
figure
plot(t,y(:,1));
xlabel('time'),ylabel('vertical displacement')
function dydt = f(t,y,hx)
k_s = 26400; %spring stiffness
m = 483; %Mass
v = 1; % speed along road
x = v*t;
hs = hfn(x,hx);
dydt =[y(2);
-( k_s*(y(1)-hs)/ m )];
end
function hs = hfn(x, hx)
hs = interp1(hx(:,1),hx(:,2),x);
end
