options = odeset('RelTol',1e-9,'AbsTol',1e-9);
[tt,yy] = ode45(@problem,[0,3],[0,0,1e-3,0]);
function z = command(t,t1,t2)
z(t>t1&t<t2)=(t(t>t1&t<t2)-t1)*t1/(t2-t1);
function dydt = problem(t,y)
data.accumulator.V_N2 = 10e-3;
data.accumulator.P_N2 = 2.5e6;
data.accumulator.p0 = 21e6;
data.accumulator.gamma = 1.2;
data.accumulator.V0 = (data.accumulator.P_N2*data.accumulator.V_N2)/(data.accumulator.p0);
data.delivery.D23 = 18e-3;
data.delivery.f23 = 0.032;
data.distributor.kd = 12;
data.distributor.d0 = 5e-3;
data.actuator.Dc = 50e-3;
data.actuator.Dr = 22e-3;
data.actuator.xmax = 200e-3;
alpha = 2*acos(1 - 2*abs(z));
A23 = (data.delivery.D23^2/4)*pi;
Ar = (data.return.Dr^2/4)*pi;
Ad = ((data.distributor.d0^2/4)*(alpha - sin(alpha)))+sqrt(eps);
A4 = (data.actuator.F0 + data.actuator.k*y(1))/ (data.accumulator.p0 - (data.tank.pT*(1-0.3)));
Peq = (data.tank.pT*A5 + (data.actuator.F0 + data.actuator.k*y(1)))/A4;
PA = data.accumulator.p0 * (data.accumulator.V0/(data.accumulator.V_N2 - y(4)))^data.accumulator.gamma;
P1 = PA - 0.5*data.delivery.ka*data.fluid.rho*((y(2)*A4)/A23)^2;
P2 = P1 - 0.5*data.delivery.ka*data.fluid.rho*((y(2)*A4)/A23)^2;
P3 = P2 - data.delivery.f23*(data.delivery.L23/data.delivery.D23)*0.5*data.fluid.rho*((y(2)*A4)/A23)^2;
P7 = data.tank.pT + 0.5*data.distributor.kd*data.fluid.rho*((y(2)*A5)/Ar)^2;
P6 = P7 + data.return.f67*(data.return.L67/data.return.Dr)*0.5*data.fluid.rho*((y(2)*A5)/Ar)^2;
P4 = P3 - 0.5*data.distributor.kd*data.fluid.rho*((y(2)*A4)/Ad)^2;
P5 = P6 + 0.5*data.distributor.kd*data.fluid.rho*((y(2)*A5)/Ad)^2;
P4 = P6 + 0.5*data.distributor.kd*data.fluid.rho*((y(2)*A5)/Ad)^2;
P5 = P3 - 0.5*data.distributor.kd*data.fluid.rho*((y(2)*A4)/Ad)^2;
dydt(2) = (1/data.actuator.ma)*((P4*A4) - (P5*A5) - (data.actuator.F0+(data.actuator.k*y(1))));
if y(1) >= data.actuator.xmax
y(1) = data.actuator.xmax;
if y(1) == data.actuator.xmax && dydt(1)>0
if (y(1) == data.actuator.xmax && dydt(2) > 0) || (y(1) == 0 && dydt(2) < 0)
if y(4) >=data.accumulator.V0