function dy=wind_4el_a(t,y,z)

% y = initial values/solved values
% z = other parameters
Il  = y(1);
AoP = y(2);
Pa  = y(3);
Pv  = y(4);
LVP = y(5);

bpm     = z(1,1);
elast   = z(:,2);
delast  = z(:,3);
etime   = z(:,4);
Rc      = z(1,5);
L       = z(1,6);
R       = z(1,7);
Ca      = z(1,8);
Cv      = z(1,9);
CAoP    = z(1,10);

if(LVP > AoP) Rao=.01;
else Rao=1000;
end

if(Pv > LVP) Rmv=.01;
else Rmv=1000;
end

te=mod(t,60/bpm); %Used to get elastance at single point in time
e=interp1(etime,elast,te); %elastance at time t
de=interp1(etime,delast,te); %de/dt at time t

Clv = 1./e;
dClv = -1./(e.^2) *de;

%     Il  AoP  Pa  Pv  LVP
a11 = [0, 1, -1, 0, 0] / L;
a21 = [-1, -(1/Rc + 1/Rao), 1/Rc, 0, 1/Rao ] / CAoP;
a31 = [1, 1/Rc, -(1/Rc + 1/R), 1/R, 0]/Ca;
a41 = [0, 0, 1/R, -(1/R + 1/Rmv), 1/Rmv]/Cv;
a51 = [0, 1/Rao, 0, 1/Rmv, -(1/Rao + 1/Rmv + dClv)]/Clv;
    
A = [a11;a21;a31;a41;a51];

dy = A*y;

