Convertir un sistema ODE de enlace simultáneo en un script de Matlab, falta de coincidencia de salida

Hola, estoy trabajando para convertir el siguiente modelo térmico en un modelo Matlab, enlace

Aquí están mis códigos para ello:

1- Función ODE

        function dX_dt = odes_thermal5(~ ,X, Tout,stat)
        %         TinIC = 26;
        r2d = 180/pi;
        Troom=X(1);
        Qlosses=X(2);
        Qheater=X(3);
        %         stat=X(4);
        % -------------------------------
        % Define the house geometry
        % -------------------------------
        % House length = 30 m
        lenHouse = 30;
        % House width = 10 m
        widHouse = 10;
        % House height = 4 m
        htHouse = 4;
        % Roof pitch = 40 deg
        pitRoof = 40/r2d;
        % Number of windows = 6
        numWindows = 6;
        % Height of windows = 1 m
        htWindows = 1;
        % Width of windows = 1 m
        widWindows = 1;
        windowArea = numWindows*htWindows*widWindows;
        wallArea = 2*lenHouse*htHouse + 2*widHouse*htHouse + ...
            2*(1/cos(pitRoof/2))*widHouse*lenHouse + ...
            tan(pitRoof)*widHouse - windowArea;
        % -------------------------------
        % Define the type of insulation used
        % -------------------------------
        % Glass wool in the walls, 0.2 m thick
        % k is in units of J/sec/m/C - convert to J/hr/m/C multiplying by 3600
        kWall = 0.038*3600;   % hour is the time unit
        LWall = .2;
        RWall = LWall/(kWall*wallArea);
        % Glass windows, 0.01 m thick
        kWindow = 0.78*3600;  % hour is the time unit
        LWindow = .01;
        RWindow = LWindow/(kWindow*windowArea);
        % -------------------------------
        % Determine the equivalent thermal resistance for the whole building
        % -------------------------------
        Req = RWall*RWindow/(RWall + RWindow);
        % c = cp of air (273 K) = 1005.4 J/kg-K
        c = 1005.4;
        % -------------------------------
        % Enter the temperature of the heated air
        % -------------------------------
        % The air exiting the heater has a constant temperature which is a heater
        % property. THeater =20 deg C
        THeater = 50;
        % Air flow rate Mdot = 1 kg/sec = 3600 kg/hr
        Mdot = 3600;  % hour is the time unit
        % -------------------------------
        % Determine total internal air mass = M
        % -------------------------------
        % Density of air at sea level = 1.2250 kg/m^3
        densAir = 1.2250;
        M = (lenHouse*widHouse*htHouse+tan(pitRoof)*widHouse*lenHouse)*densAir;
        %%
        %%
        dQheater_dt= (THeater-Troom)*Mdot*c*stat;
        dQlosses_dt=(Troom-Tout)/Req;
        dTroom_dt=(1/(M*c))*( dQheater_dt - dQlosses_dt );
        dX_dt=[dTroom_dt;dQlosses_dt;dQheater_dt];
        end

Y así es como lo llamé:

    % clear;
    % t=0:1:48;
    % f=1/3600;
    % Tout = @(t) T_const + k * sin(t);
    % % Tout =  k * sin(2*pi*f*t);
    % plot(t,Tout(1:end))
    clear;
    T_const = 50; k = 15;
    f=2/48;
    ts=1/3600;
    T=48;
    tsin=1:ts:T;
    y=T_const+k*sin(2*pi*f*tsin);
    Tout = decimate( y , 3599 );
    tsin2 = decimate( tsin, 3599 );
    t1=22;t2=26;
    Troom=zeros(48,1);
    Qheater=zeros(48,1);
    Qlosses=zeros(48,1);
    stat=zeros(48,1);
    %
    Troom_vec=[];
    Qheater_vec=[];
    Qlosses_vec=[];
    tm=[];
    for i=1:48
        if(i==1)
            Troom_0=28;Qlosses_0=0;Qheater_0=0;
            X0=[Troom_0;Qlosses_0;Qheater_0];
            stat(i)=1;
        end
        if(i==2)
            Troom_0=Troom(i-1);Qlosses_0=Qlosses(i-1);Qheater_0=Qheater(i-1);
            X0=[Troom_0;Qlosses_0;Qheater_0];
            stat(i)=1;
        end
        if(i>2)
            Troom_0=Troom(i-1);Qlosses_0=Qlosses(i-1);Qheater_0=Qheater(i-1);
            X0=[Troom_0;Qlosses_0;Qheater_0];
            if(Troom(i-1)-Troom(i-2)>0 & Troom(i-1)>t1 & Troom(i-1)<t2)
                stat(i)=1;
            end
            if(Troom(i-1)-Troom(i-2)<=0 & Troom(i-1)>t1 & Troom(i-1)<t2)
                stat(i)=0;
            end
            if(Troom(i-1)>t2)
                stat(i)=0;
            end
            if(Troom(i-1)<t1)
                stat(i)=1;
            end
        end
        tspan = i-1:0.1:i;  % Obtain solution at specific times
        %      tspan = 0:1;  % Obtain solution at specific times
        [tx,X] = ode45(@(t,y) odes_thermal5(t,y,Tout(i),stat(i)),tspan,X0);
        Troom_vec=[Troom_vec;X(:,1)];
        Qlosses_vec=[Qlosses_vec; X(:,2)];
        Qheater_vec= [Qheater_vec;X(:,3)];
        Troom(i)=X(6,1);
        Qlosses(i)=X(6,2);
        Qheater(i)=X(6,3);
        tm=[tm;tx];
    end
    %         stat=X(:,4);
    figure(1);
    plot(tsin2,Troom);
    hold on
    plot(tsin2,stat*10);
    hold off
    ylabel('Troom');
    xlabel('t');

Mi salida es básicamente Troom, pero la respuesta que obtuve fue Tout, como si no tuviera un modelo de calentador. Aquí es cómo se ve la salida, ¿alguna idea de por qué es así?