% 5-element Winkessel Model

V1 = 12.2;
I2 = 0.0078;
V3 = 12.0;

C1 =   0.002;
L2 =   0.05;
C3 =   0.015;
Rd1 = 1e6;
R2 = 4;
R3  = 220;

Data = zeros(3000,5);
t0 = zeros(3000,1);

t = 0;
dt = 1e-5;
npt = 0.001/dt;

for i=1:3000

   for j=1:npt

      V0 = max(0, 17.6*sin(6.28*t));

      if (V0 > V1)  Rd1 = 0.1;  
               else Rd1 = 1e6;
               end
      
      I0 = (V0 - V1)/Rd1;

      dV1 = ( I0 - I2 )/C1;
      dI2 = (V1 - R2*I2 - V3)/L2;
      dV3 = (I2 - V3/R3)/C3;

      t = t + dt;
      V1 = V1 + dV1*dt;
      I2 = I2 + dI2*dt;
      V3 = V3 + dV3*dt;

      end
    
    t0(i) = t;
    Data(i,:) = [I0*1000, V0, V1, I2*1000, V3];

   end

clf
subplot(211);

plot(t0, Data(:,[2,3,5])*7.5);
xlabel('Time (seconds)');
ylabel('Pressure (mmHg)');
xlim([1,3]);

subplot(212);
plot(t0,Data(:,[1,4]))
xlabel('Time (seconds)');
ylabel('Flow (ml/s)');
xlim([1,3]);