Simulink thermal model in Matlab

Hello I am working on converting this model into a Matlab model, Here are my codes for it:

1- ODE function

            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

And here is how I called it:

        % 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');

My output is basically Troom, but the response I got was Tout, as if I don't have a heater model. Here is how the output looks like, any idea why its like that?

Thanks