MATLAB: Convert transfer function to state-space models

Control System ToolboxlsimMATLABstate-spacetransfer function

I'm trying to convert a transfer function to state space model. In order to avoid the control toolbox functions since I don't have it's license.
The transfer function that I'm using is quite simple. Is the model for a heatsink temperature.
For this transfer function I have the next values for A,B,C and D.
I've tried to compare the output from the lsim() function and the equations from above but the output differs
Here's the code that I've used
P = rand(1,10)*1000;
t = 0:1:length(P)-1;
A = -400;
B = 0.5;
C = 0.6658;
D = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%



%%%    State space    %%%
%%%%%%%%%%%%%%%%%%%%%%%%%
state_space_sys = ss(A,B,C,D);
state_space = lsim(state_space_sys,P,t);
figure
plot(t,state_space)
grid on,grid minor, title('State space')
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   My function     %%%
%%%%%%%%%%%%%%%%%%%%%%%%%
x(1:length(P)) = 0;
y(1:length(P)) = 0;
u = P;
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(t,y)
grid on,grid minor, title('My function')
Is there any mistakes on the approach?
It should be quite simple, but I can't manage to find the solution.

Best Answer

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
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() and the function lsim() with the A,B,C and D values.
Related Question