tspan = linspace(0, 1.2, 12001);
[t, x] = ode45(@clutchStateFcn, tspan, x0);
[~, u(j)] = clutchStateFcn(t(j), x(j,:).');
r_slip = max(0, 70 - 70*(t + 0.45));
v_slip = x(:,1) - x(:,2);
tL = tiledlayout(2, 1, 'TileSpacing', 'Compact');
plot(t, u), grid on, ylim([-1000, 7000])
title('Open-loop Discontinuous Clamping Force (Input)')
plot(t, [r_slip, v_slip]), grid on
legend('Reference', 'Actual')
xlabel(tL, 'Time / s'), title('Slip Velocity (Output)')
text(0.10, 60, '\leftarrow 0.0 kN')
text(0.32, 40, '\downarrow 1.5 kN')
text(0.55, 20, '\downarrow 4.0 kN')
text(0.93, 20, '\downarrow 6.0 kN')
function [dxdt, u] = clutchStateFcn(t, x)
ds = 2*zeta*sqrt(ks*jd*jv/(jd+jv));
g1 = (np*Rm*(mu0 + muk*Vslip)*tanh(2*Vslip))/je;
g2 = (np*Rm*(mu0 + muk*Vslip)*tanh(2*Vslip))/jd;
f2 = - ds/jd*x(2) + ds/jd*x(3) + ks/jd*x(4);
r_Fnn(t >= 0.3 & t < 0.45) = 1500;
r_Fnn(t >= 0.45 & t < 0.8) = 4000;