% Define muscle parameters
m = 1;  % mass of the muscle
b = 10;  % damping coefficient
k = 100;  % stiffness coefficient
l0 = 1;  % initial muscle length

% Define stretch reflex parameters
T = 0.1;  % delay time of the reflex loop
G = 10;  % gain of the reflex
gamma = 1;  % gamma motoneuron gain

% Define input stretch
t = 0:0.01:1;  % time vector
x = 0.1 * ones(size(t));  % constant stretch of 0.1 units

% Initialize variables
y = zeros(size(t));  % muscle length
v = zeros(size(t));  % muscle velocity

% Initial conditions
y0 = l0;
v0 = 0;

% Loop through time steps
for i = 2:length(t)
    % Handle edge cases for delay
    if i <= T / 0.01 + 1
        reflex_activation = gamma * (l0 - y(i-1));
    else
        reflex_activation = G * (x(i-round(T/0.01)) - y(i-round(T/0.01))) + gamma * (l0 - y(i-1));
    end

    % Muscle force calculation
    f = k * (y(i-1) - x(i-1)) - b * v(i-1);

    % Muscle acceleration
    a = (f - reflex_activation) / m;

    % Update muscle velocity and length
    v(i) = v(i-1) + a * (t(i) - t(i-1));
    y(i) = y0 + v(i) * (t(i) - t(i-1));
end

% Plot results
plot(t, y, 'b');
xlabel('Time (s)');
ylabel('Muscle Length (units)');
title('Steady State Muscle Length Response');