%%
% Problem constants
m = 1;
c = 2;
k = 5;
g = 9.8;

% Time vector for ODE45
tstep = 0.01;
tspan = 0:tstep:20;

% seismic
%---with function
%groundfunc = @(t) ground_step_tanh(t);   %Seismic
%ground = groundfunc(tspan);
%---random ground
%ground = smooth(cumsum(normrnd(0,0.005,length(tspan),1)),10);
%---flat ground, y(t) = 0
%ground = zeros(length(tspan),1);
%---explicit
ground = 0.1*sin(5*t);
% y prime
groundp = gradient(ground,t);

% ODE function
odefunc = @(t, x) bot_susp_1d(t, x, m, k, c, tspan, ground, groundp); %for ode45 syntax

% State vector: [x(t) x'(t)]
initial_x = - m * g / k;
initial_xp = 0;

% initial_yp = 0;
initial = [initial_x initial_xp];

%% ode45 execution
tic
[t,xres] = ode45(odefunc,tspan, initial);
toc

%% Plotting
plot(t, xres(:,1));