function ode = model(x,u,d,p)
V = x(1);                                       % [m3] Hold up CSTR reactor
ca = x(2);                         % [kmol/m3] Concentration of component A

q1 = u(1);                                          % [m3/min] Inlet stream 
q2 = u(2);                                         % [m3/min] Outlet stream 

caf = d(1);                % [kmol/m3] Concentration of component A in feed
k = d(2);                            % [m3/kmol/min] Reaction rate constant  

%% model equations

dVdt  = q1-q2;

dcadt = (q1/V)*(caf-ca)-k*ca^2;

error = 0;
if (dVdt == inf || dcadt == inf)
    dVdt = 0;
    dcadt =0;
    error =-1;
end
if (isnan(dVdt) || isnan(dcadt))
    dVdt = 0;
    dcadt =0;
    error =-1;   
end

ode = [dVdt;dcadt;error];                                           % model

end
