How to avoid tf() command using State-space models in Matlab

I'm trying to avoid the function tf() from Matlab since it requires specific toolboxes to run.

The transfer function that I'm using is quite simple. Is the model for a heatsink temperature.

                         H(s) = (Rth/Tau)/(s + 1/Tau)

In order to avoid the tf() function, I've tried to substitute the transfer function with a state space model coded in Matlab.

I've used the function ss() to get te values of A,B,C and D. And I've tried to compare the results from tf() and my function.

Here's the code that I've used:

Rth = 8.3220e-04; % ºC/W
Tau = 0.0025; % s

P = rand(1,10)*1000; % Losses = input
t = 0:1:length(P)-1; % Time array

%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer function %%%
%%%%%%%%%%%%%%%%%%%%%%%%%

H = tf([0 Rth/Tau],[1 1/Tau]);
Transfer_func = lsim(H,P,t);

figure, plot(Transfer_func),grid on,grid minor, title('Transfer func')

%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   My función ss   %%%
%%%%%%%%%%%%%%%%%%%%%%%%%

% Preallocate for speed
x(1:length(P)) = 0;
y(1:length(P)) = 0;
u = P;

sys = ss(H);
A = sys.A;
B = sys.B;
C = sys.C;
D = sys.D;

for k = 1:length(u)
     x(k+1) = A*x(k) + B*u(k);
     y(k) = C*x(k) + D*u(k);
end

figure, plot(y), grid on,grid minor, title('With my función')

I know that the values from A,B,C and D are ok, since I've checked them using

H = tf([0 Rth/Tau],[1 1/Tau]);
sys = ss(H);
state_space_sys = ss(sys.A,sys.B,sys.C,sys.D);
state_space = lsim(state_space_sys,P,t);

figure, plot(state_space),grid on,grid minor, title('State space')

As you can see, the results obtained from my funtion and the function tf() are very different.

Is there any mistakes on the approach?

If it's not possible to avoid the tf() function in this way, is there any other way?

Solution

At the end, I found another solution. I'm posting this here, so if someone has the same problem, can use this approach.

If you take the transfer function, and develop it, we reach to the following expresion

H(s) = deltaT(s)/P(s) = (Rth/Tau)/(s + 1/Tau)

deltaT(s) * (s + 1/Tau) = (Rth/Tau) * P(s)

deltaT(s) * s = (Rth/Tau) * P(s) - deltaT(s)/Tau

Now, we know that 1/s is equal to integrate. So in the end, we have to integrate the right side of the equation. The code would be like this.

Cth = Tau/Rth;
deltaT = zeros(size(P));

for i = 2:length(P)
    deltaT(i) = (1/Cth * (P(i)-deltaT(i-1)/Rth))*(time(i)-time(i-1)) + deltaT(i-1);
end

This integral has the same output as the function tf().