clear; close all
%
%  Spring-damper system, initial displacement
%
k = 10000;  %% N/m
b = 1000;   %% N-sec/m
tc = k/b;   %% time constant
x0 = 0.1;   %% initial displacement
duration = 5*tc;
dt = duration/200;
t = 0:dt:duration;
%
%  Solve for y(t) and y'(t)
%
%  bksys = name of our function to get y' given y,t
%  [0 5*tc] = start & end times
%  y0 = initial y and y'
%  b & k passed to mbksys by ode45
%
[t,x] = ode45(@bksys,t,x0,[],b,k);
%
%  Plot x(t)
%
subplot(2,1,1);
plot(t,x(:,1)); grid
xlabel('Time');
ylabel('Displacement')
title('\bf spring-damper system with initial displacement')
%
%  Get x'(t), plot
%
xp = bksys(t,x,k,b);
subplot(2,1,2);
plot(t,xp); grid
xlabel('Time');
ylabel('Velocity');
%
%  Plot forces
%
figure
plot(t,k*x,t,b*xp);
legend('Spring force','Damper force'); grid