ODE45 outputs only zeros

I'm trying to use ode45 to model the movement of a piston. This works though a command (function command) that opens with respect of time. When I try to run the code the only output I get is zeros even if the output I expect for yy(:,1) should be something like the image attached.

I don't understand what's wrong with my code and why it outputs zeros instead of numbers.

t0 = 0;
t1 = 1;
t2 = 1.5;
tf = 3;

options = odeset('RelTol',1e-9,'AbsTol',1e-9);
[tt,yy] = ode45(@problem,[0,3],[0,0,1e-3,0]);
plot(tt,yy(:,1))

function z = command(t,t1,t2) %FUNCTION FOR COMMAND Z
   z=zeros(size(t));
   z(t>t1&t<t2)=(t(t>t1&t<t2)-t1)*t1/(t2-t1);
   z(t>=t2)=t1;
end


function dydt = problem(t,y)

%------------------------------------------------------------------------------------%
data.fluid.rho = 890;

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.ka = 1.12;
data.delivery.kcv = 2;
data.delivery.D23 = 18e-3;
data.delivery.L23 = 2;
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.ma = 2;
data.actuator.xmax = 200e-3;
data.actuator.F0 = 1e3;
data.actuator.k = 120e3;

data.return.Dr = 18e-3;
data.return.L67 = 15;
data.return.f67 = 0.035;

data.tank.pT = 0.1e6;
data.tank.VT = 1e-3;
data.tank.kt = 1.12;
%------------------------------------------------------------------------------------%
t1 = 1;
t2 = 1.5;

z = command(t,t1,t2);
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)));
A5 = (1-0.3)*A4;
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 = PA - (data.fluid.rho*(y(2)*A4)^2)*0.5*( (data.delivery.ka/A23^2) + (data.delivery.kcv/A23^2) + (data.distributor.kd/Ad^2) + (data.delivery.f23*(data.delivery.L23/(data.delivery.D23*A23^2))));
% P5 = data.tank.pT + (data.fluid.rho*(y(2)*A5)^2)*0.5*( (data.tank.kt/Ar^2) + (data.distributor.kd/Ad^2) +   (data.return.f67*(data.return.L67/(data.return.Dr*Ar^2))));


if z == 0
    P4 = Peq;
    P5 = data.tank.pT;
elseif z > 0
    if Ad < 1e-9
        P4 = Peq;
        P5 = data.tank.pT;
    else
        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;
    end
else
    if Ad < 1e-9
        P4 = Peq;
        P5 = data.tank.pT;
    else
        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;
    end
end

dydt = zeros(4,1);
%y(1) position; y(2) velocity
dydt(1) = y(2);  %  xa_dot = va; 

dydt(2) = (1/data.actuator.ma)*((P4*A4) - (P5*A5) - (data.actuator.F0+(data.actuator.k*y(1))));  %  va_dot = (1/ma)*((P4*A4) - (P5*A5) - (F_0+(k*xa)));

dydt(3) =  y(2)*A5;  %  Vt_dot = Qout;
 
dydt(4) =  y(2)*A4;  %  Vacc_dot = -Qin;


%BOUNDARY CONDITION FOR THE PISTON
% position cant go below zero nor xmax
%if velocity is positive when x=xmax put to zero
%if negative is positive when x=xmax put to zero
%acceleration must go to zero when piston arrives at x=0 and x=xmax
if y(1) <=0
    y(1)=0;
end
if y(1)==0 && dydt(1)<0
    dydt(1)=0;
end
if y(1) >= data.actuator.xmax
    y(1) = data.actuator.xmax;
end
if y(1) == data.actuator.xmax && dydt(1)>0
    dydt(1) = 0;
end



if (y(1) == data.actuator.xmax && dydt(2) > 0) || (y(1) == 0 && dydt(2) < 0)
    dydt(1) = 0;
    dydt(2) = 0;
end

% when the volume of liquid inside the accumulator is bigger than the
% initial one put to zero
if y(4) >=data.accumulator.V0
    dydt(4) = 0;
end








end

Answer 1

[Bjorn Gustavsson]

When you have discontinuities in the ODEs (the derivative of command) then the safe approach is to separate the integration into sections where everything is continuous. You might get somewhere by reducing the max-stepsize of the ode-solever by setting MaxStep in the odeoptions-struct returned from odeset to something smaller, but really think about integrating in chunks where everyting behaves nicely.

HTH