function out = model(x,u,d,p)
% model for level and concentration 
T      = x;                                  % [K] temperature of cold side

z      = u;                                % [-] steam inlet valve position

T0     = d(1);                            % [K] cold side inlet temperature
w      = d(2);                                      % [kg/s] cold side flow
pd0    = d(3);                            % [Bar] steam uppstreams pressure

m      = p(1);                                  % [kg] hold up og cold side
cp     = p(2);                        % [kJ/kgK] heat capacity of cold side
U      = p(3);                         % [kW/m2K] heat transfer coefficient 
A      = p(4);                                 % [m2] Area of heat transfer
lambda = p(5); % [kJ/kg] The enthalpy of vaporization of medium in hot side
k      = p(6);                               % [kg/Bar*s] valve coefficient
R      = p(7);                              % [kJ/kmolK] ideal gas constant
pdref  = p(8);                  % [Bar] saturated steam refference pressure
Tdref  = p(9);                 % [K] saturated steam refference temperature
Mm     = p(10);                             % [Kg/kmol] Molar mass of steam
% model equation

% need to solve a implicite system of equations to find Td (use a help
% function)
Rs = R/Mm;

pd = @(Tdvar) pdref*exp(-(lambda/Rs)*(1/Tdvar -1/Tdref));

fun = @(Tdvar) k*z*lambda/(U*A)*(pd0-pd(Tdvar)) + T - Tdvar;

options = optimset('Display','off','TolFun',1e-20,'TolX',1e-20,'MaxIter',1e4,'MaxFunEvals',1e10);

Td  = fsolve(fun,Tdref,options);

Q = U*A*(Td - T);

wd = Q/lambda;

dTdt = (1/m)*(w*(T0-T)+Q/cp);

out = [dTdt,wd,Td];                                                 % model

end