why do i get a wrong answers when using the if else sentences in the defined-function?

hi all! I use a new method to solve ode numerically, called QSS, it doesn't matter if you can't understand the method. My problem is that when I vectorized my codes which means I have to defined the ode in the user defined-function, bacause the ode is nonsmooth, I have to use the 'if else' sentence in the defined-function, however the results is not right compared to the original codes(the original codes puts the if else sentence in the loop), can anyone tell me how could that happen? Why dosen't the if else sentence that put in the function work? I really need someone to help me fix the problem. Thank you all:)

Below is my codes.First one is the original codes, second one is the vectorized codes.

Besides, both codes work well when solving the smooth odes.

clc
clear
tic
%为各变量赋初值
delta_q=1e-4;
q1=21;q2=25;q3=17;
x1=q1;x2=q2;x3=q3;
t=0;delta_t=0;
A=zeros(105609,6);
n=0;
C1=550;C2=600;Cw=580;R1=1.8;R2=2.2;R1w=2;R2w=2;Ta=32;p1=6;p2=6;m1=1;m2=1;T_ref=20;pi=3.1415926;
%开始while循环
while (t<3600*12)
    if(x1>T_ref+0.5)
        m1=1;
    elseif(x1<T_ref-0.5)
        m1=0;
    end
    if(x2>T_ref+0.5)
        m2=1;
    elseif(x2<T_ref-0.5)
        m2=0;
    end
     Ta=10*sin(pi*3.4722e-05*t-pi/2)+20;
     TTa=pi*3.4722e-05*10*cos(pi*3.4722e-05*t-pi/2);
    TTTa=-(pi*3.4722e-05)^2*10*sin(pi*3.4722e-05*t-pi/2);
    Dx1=((Ta-q1)/R1+(q3-q1)/R1w-m1*p1)/C1;  
    Dx2=((Ta-q2)/R2+(q3-q2)/R2w-m2*p2)/C2; 
    Dx3=((q1-q3)/R1w+(q2-q3)/R2w)/Cw;
    DDx1=((TTa-Dx1)/R1+(Dx3-Dx1)/R1w)/C1;
    DDx2=((TTa-Dx2)/R2+(Dx3-Dx2)/R2w)/C2;
    DDx3=((Dx1-Dx3)/R1w+(Dx2-Dx3)/R2w);
   %求Δt1和Δt2的值
       delta_t1=sqrt(2*delta_q/abs(DDx1));
       delta_t2=sqrt(2*delta_q/abs(DDx2));
       delta_t3=sqrt(2*delta_q/abs(DDx3));
       
     delta_tmin=min([delta_t1 delta_t2 delta_t3]);


   %比较Δt1和Δt2的大小，进而确定q1和q2谁先跃迁
   if (delta_t1==delta_tmin)
       delta_t=delta_t1;
       t=t+delta_t;
       x1=x1+Dx1*delta_t+0.5*DDx1*delta_t^2;
       x2=x2+Dx2*delta_t+0.5*DDx2*delta_t^2;
       x3=x3+Dx3*delta_t+0.5*DDx3*delta_t^2;
       q1=x1;
       q2=q2+Dx2*delta_t;
       q3=q3+Dx3*delta_t;
   elseif (delta_t2==delta_tmin)
       delta_t=delta_t2;
       t=t+delta_t;
       x1=x1+Dx1*delta_t+0.5*DDx1*delta_t^2;
       x2=x2+Dx2*delta_t+0.5*DDx2*delta_t^2;
       x3=x3+Dx3*delta_t+0.5*DDx3*delta_t^2;
       q2=x2;
       q1=q1+Dx1*delta_t;
       q3=q3+Dx3*delta_t;
   elseif (delta_t3==delta_tmin)
       delta_t=delta_t3;
       t=t+delta_t;
       x1=x1+Dx1*delta_t+0.5*DDx1*delta_t^2;
       x2=x2+Dx2*delta_t+0.5*DDx2*delta_t^2;
       x3=x3+Dx3*delta_t+0.5*DDx3*delta_t^2;
       q3=x3;
       q1=q1+Dx1*delta_t;
       q2=q2+Dx2*delta_t;
   end
   n=n+1;
   A(n,1)=t;
   A(n,2)=x1;
   A(n,3)=x2;
   A(n,4)=x3;
   A(n,5)=m1;
   A(n,6)=m2;
end
tt=0:3600*12;
yy=10*sin(pi*3.4722e-05*tt-pi/2)+20;
figure
  plot(A(:,1),A(:,2),'--','LineWidth',1.25);hold on
  plot(A(:,1),A(:,3),'-.','LineWidth',1.25);
  plot(A(:,1),A(:,4),'-','LineWidth', 1.25);
  plot(tt,yy,'-k');
  xticks(0:3600:13*3600); %这样x轴会每隔10显示一个刻度
  xticklabels({'6','7','8','9','10','11','12','13','14','15','16','17','18'});%为说明效果，省略了部分内容，写代码时应该与xticks对应写够10+1个label
  xlabel('时刻（h）','FontSize',12,'FontWeight','bold','Color','k');
  ylabel('摄氏度（℃）','FontSize',12,'FontWeight','bold','Color','k');
  grid on
   legend('room1','room2','wall','ambient temperature')
% figure
%    plot(A(:,1),A(:,5));hold on
%    plot(A(:,1),A(:,6));
toc     

second one

% QSS1矢量化
clc
clear
tic
%为各变量赋初值
delta_q=[1; 1; 1]*1e-4;
q=[21; 25; 17];
x=q;
t=0;
A=zeros(107666,length(q)+1);
delta_x=zeros(length(q),1);
n=0;
%开始while循环
while (t<3600*12)
    
    ff=Dx(t,q);
    aa=ff(1:length(q));
    aaa=ff(length(q)+1:end);
   %求Δt1和Δt2的值
    delta_t=sqrt(2*delta_q./abs(aaa));
   %比较Δt1和Δt2的大小，进而确定q1和q2谁先跃迁
    cc=find(delta_t==min(delta_t));
   %存在多个最小值相等，取第一个最小值
    if length(cc)>1
        cc=cc(1);
    end
    dd=find(delta_t~=min(delta_t));
    t=t+delta_t(cc);
    x=x+aa.*delta_t(cc)+0.5*aaa.*delta_t(cc)^2;
    q(cc)=x(cc);
    q(dd)=q(dd)+aa(dd).*delta_t(cc);
   %赋值 循环
    n=n+1;
    A(n,1)=t;
    A(n,2)=x(1);
    A(n,3)=x(2);
    A(n,4)=x(3);
end
toc

% % 绘图
hold all;
LineWidth=1.25;
tt=0:3600*12;
yy=10*sin(pi*3.4722e-05*tt-pi/2)+20;
plot(A(:,1),A(:,2),'--','LineWidth',LineWidth);
plot(A(:,1),A(:,3),'-.','LineWidth',LineWidth);
plot(A(:,1),A(:,4),'-','LineWidth', LineWidth);
plot(tt,yy,'-k');
xticks(0:3600:13*3600); %这样x轴会每隔10显示一个刻度
xticklabels({'6','7','8','9','10','11','12','13','14','15','16','17','18'});%为说明效果，省略了部分内容，写代码时应该与xticks对应写够10+1个label
  xlabel('时刻（h）','FontSize',12,'FontWeight','bold','Color','k');
  ylabel('摄氏度（℃）','FontSize',12,'FontWeight','bold','Color','k');
grid on; grid minor;
% title({'QSS量子=',string(delta_q(1)) });
% legend
% legend('QSS求解');
%比较误差

function f=Dx(t,x)
    C1=550;C2=600;Cw=580;R1=1.8;R2=2.2;R1w=2;R2w=2;
    p1=6;p2=6;m1=1;m2=1;T_ref=20;pi=3.1415926;
    if(x(1)>T_ref+0.5)
        m1=1;
    elseif(x(1)<T_ref-0.5)
        m1=0;
    end
    if(x(2)>T_ref+0.5)
        m2=1;
    elseif(x(2)<T_ref-0.5)
        m2=0;
    end
     Ta=10*sin(pi*3.4722e-05*t-pi/2)+20;
     TTa=pi*3.4722e-05*10*cos(pi*3.4722e-05*t-pi/2);
    f=zeros(6,1);
    f(1)=((Ta-x(1))/R1+(x(3)-x(1))/R1w-m1*p1)/C1;
    f(2)=((Ta-x(2))/R2+(x(3)-x(2))/R2w-m2*p2)/C2; 
    f(3)=((x(1)-x(3))/R1w+(x(2)-x(3))/R2w)/Cw;
    f(4)=((TTa-f(1))/R1+(f(3)-f(1))/R1w)/C1;
    f(5)=((TTa-f(2))/R2+(f(3)-f(2))/R2w)/C2;
    f(6)=((f(1)-f(3))/R1w+(f(2)-f(3))/R2w);
end

[Walter Roberson]

http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F

Do not compare calculated numbers for equality.

Also, even at the best of times, all three of those if/elseif tests are going to fail. I suspect that you are looking to test whether the values are in particular ranges, which requires a range test rather than an equality test.

if value <= first_boundary
    do first range stuff
elseif value <= second_boundary
    do second boundary stuff
elseif value <= third_boundary
    do third boundary stuff
else
    do stuff when you are further than any bound
end

[汉武 沈]

Thank you for your comment.I am still confused.The first above codes is the original one, which compare for equilty works fine, the codes is okay. The second one failed, it doesn't compare for equilty. My problem in the second codes is that I put the 'if else' in the defiined-function which will be passed to the main function just like we did when we use the ode suits, however it seemed that the 'if else' sentence didn't work, the variables x(1)、x(2) should fluctuate between 20.5 and 19.5, if you run the first codes, you will get the right answers.Is there something wrong when 'if else' sentences be put in the defined-function? As you can see the first codes also catains the exactly same 'if else' sentences, but it works well.