function [sys] = cortinez(Age,Weight,stepSize)
V1std = 4.48;
V2std = 21.2;
V3std = 237;
CL1std = 1.92;
Q2std = 1.45;
Q3std = 0.86;
SLV2 = -0.0164;
SLQ2 = -0.0153;
wa = Weight/70;
wap = wa^0.75;
aa = Age - 50;
V1 = V1std * wa;
V2 = V2std * wa * exp(SLV2*aa);
V3 = V3std * wa;
volume = [V1 V2 V3];
CL1 = CL1std* wap;
CL2 = Q2std* wap * exp(SLQ2*aa);
CL3 = Q3std*wap;
clearance = [CL1 CL2 CL3];
tpeak = 1.6;
ike0 = [.1;2];
t=0:.0005:3;
u=zeros(size(t));
u(1)=1;
options = optimset('Display','off');
ke0 = fzero(@getKe0, ike0, options);
function [error] = getKe0(tke0)
[sys] = mam2ss(clearance,volume,tke0);
[y,t] = lsim(sys,u,t);
s = lsiminfo(y,t,0);
error = tpeak - s.MaxTime;
end
sys = mam2ss(clearance,volume, ke0);
if (nargin==3)
sys=c2d(sys, stepSize,'zoh');
end
end
function [sys] = mam2ss(clearance,volume, ke0)
n = length(clearance);
if (nargin == 2)
k1 = clearance./volume(1);
k2 = clearance(2:n)./volume(2:n);
A = [-sum(k1) k2;transpose(k1(2:n)) -diag(k2)];
B = [1;zeros(n-1,1)];
C = [1/volume(1) zeros(1,n-1)];
else
EFFECT_VOL_FACTOR=10000;
k1 = [clearance./volume(1) ke0/EFFECT_VOL_FACTOR];
k2 = [clearance(2:n)./volume(2:n) ke0];
A = [-sum(k1) k2;transpose(k1(2:n+1)) -diag(k2)];
B = [1;zeros(n,1)];
C = [zeros(1,n) EFFECT_VOL_FACTOR/volume(1)];
end
C = C./1000;
D=0;
sys = ss(A,B,C,D);
end