% syms t;
m1=32;m2=16;c1=1013.4;K1=33000;K2=185000;f1=1000;f2=-1000;f3=500;
A=[0 1 0 0;-K1/m1 -c1/m1 K1/m1 c1/m1;0 0 0 1;K1/m2 c1/m2 -(K1+K2)/m2 -c1/m2];
B=[0;-1/m1;0;1/m2];
E=[0;0;0;K2/m2];
C=[1 0 0 0];
H=[f1 f3 f2 0];
CONZ=[0 1 0 0];
CONV=[0 1 0 -1];
CONS=[1 0 -1 0];
x=[0;0;0;0];
T=0.001;
hh=T;
N=10000;
for n = 1: N
fin(1,n) = 181/99900*randn(1,1)-181/99900*randn(1,1)*exp(-999/100*(n-1)*T);
f=fin(1,n);
CONT1 = CONV* x;
CONT2= CONZ* x;
if CONT1.* CONT2 >= 0
AA= A;
else
AA= A + B* H;
end
k0=AA* x + E* f;
k1=AA* (x+ hh* k0 / 2) + E* f;
k2=AA* (x+ hh* k1 / 2) + E* f;
k3=AA* (x+ hh* k2) + E* f;
x=x +(k0 +2* k1+ 2* k2 +k3)* hh/6;
y=C* x;
a1=-K1/m1*CONS* x-c1/m1*CONV* x;
yy(1,n) = a1;
end
t=[1:N]* T;
figure(1);
plot(t,yy);