Bs = 6000;
mc = 1500 / 4;
mw = 25;
Ks = 50000;
Kt = 90000;

%

A = [-Bs / mc, Bs / mc, 1 / mc, 0;
    Bs/ mw, -Bs / mw, -1 / mw, 1 / mw;
    -Ks, Ks, 0, 0;
    0, -Kt, 0, 0];
B = [0; 0; 0; Kt];
C = [1 0 0 0;
    0 1 0 0];
D = [0; 0];

%

dt = 0.001;
tf = 1.5;
t = 0:dt:tf;
X = zeros(4, length(t));
Y = zeros(2, length(t));
U = 0.01;

%

for i = 1:length(t)-1
    X(:,i+1) = X(:,i) + dt * (A * X(:,i) + B * U);
    Y(:,i+1) = C * X(:,i+1) + D * U;
end

%

figure()
plot(t, Y(1,:))
hold on
plot(t, Y(2,:))
legend("Car", "Wheel")
ylabel("Velocity")
xlabel("Time")
hold off

%

sys = ss(A, B, C, D);

%

opts = stepDataOptions("StepAmplitude", U);
step(sys, opts)