function Tr_sim = simulate(par)

global t u_ref y_ref data

R_h=par(1);
R_ra=par(2);
C_r=par(3);
C_h=par(4);
A_w=par(5);

% R_h=10;
% R_ra=10;
% C_r=1000;
% C_h=1000;
% A_w=2;

% Dynamikmatrix / Zustandsmatrix A von Modell 2R2C
A = [((-R_h-R_ra)/(R_ra*R_h*C_r)), 1/(R_h*C_r);...
    1/(R_h*C_h), (-1)/(R_h*C_h)];
% Eingangsmatrix B von Modell 2R2C
B = [ 1/(R_ra*C_r), A_w/C_r, 0; 0, 0, 1/C_h ];
% Ausgangsmatrix C
C = [1, 0];
% Durchgangsmatrix D
D = zeros(1,3);

sys_ss=ss(A,B,C,D); % State Space Modell

Tr_sim = zeros(length(t),1);
Tr_sim(1,1) = y_ref(1);
Tr0_sim = Tr_sim(1,1);

for i = 1:length(t)-1
    ts = [t(i),t(i+1)];
    opts = odeset('Mass',@M,'MStateDependence','none');
    sol = ode15s(@(t,u_ref)sys_ss(t,u_ref,par),ts,Tr0_sim,opts);
    Tr0_sim = sol.y(1,end);
    Tr_sim(i+1,1) = Tr0_sim(1);
end

end