%
% �¶‘ÌƒVƒXƒeƒ€�HŠw‚ÌŠî‘b
% �}4.21
%
clear all
close all
L0d = 3.57;
l0c = 0.97;
Kd = 1.24;
Kc = 0.047;
D = 3.78;
fp = 32.25;
fs = 9.43;
P0 = -1.007;
taup = 0.22;
taus = 0.77;

r(1) = 2.62;
v(1) = 0;
step = 5001;
dt = 0.001;
for i=2:step
    if (i > 1000 + (taup / dt) && i <= 1100 + (taup / dt))
        f1 = fp;
    else
        f1 = 0;
    end
    if (i > 1000 + (taus / dt) && i <= 1100 + (taus / dt))
        f2 = fs;
    else
        f2 = 0;
    end
    v(i) = v(i-1) + dt * (-Kc * (l0c - r(i-1))^2 + Kd*(L0d - r(i-1))^2 - D * v(i-1) - f1 + f2 + P0);
    r(i) =  r(i-1) + dt * v(i);
end
plot(r);
