function [C Ceq] = twocol_dyn_constr(u)

global d x_steady Vmax Vmax_2 config

%Simulate for 5000 seconds

[tsol,xsol] = ode15s(@twocol_dyn,[0 5000],x_steady,[],u,d);

%Take x vector at final time step and retrieve necessary values from model:

%[xprime D1 B1 L1 V1 xD1 xB1 e1 D2 B2 L2 V2 xD2 xB2 e2] = twocol_dyn(xsol(end,:),u,d) ;
[xprime D1 B1 L1 V1 xD1 xB1 e1 D2 B2 L2 V2 xD2 xB2] = twocol_dyn(xsol(end,:),u,d) ; %Don't need e2 for now.
if max(xprime) > 1e-6
    [tsol,xsol] = ode15s(@twocol_dyn,[5000 25000],xsol(end,:),[],u,d);
    [xprime D1 B1 L1 V1 xD1 xB1 e1 D2 B2 L2 V2 xD2 xB2] = twocol_dyn(xsol(end,:),u,d) ; %Don't need e2 for now.
    if max(xprime) > 1e-6
    disp('Steady state not reached after 25000 seconds')
    end
end

switch config
    case('XdXb')
        Cmin = [-D1 ; -B1 ; -L1 ; -V1 ; -D2 ; -B2 ; -L2 ; -V2] ; %Make sure all flows are positive!
        Cmax = [V1 - Vmax ; V2 - Vmax_2] ;
        C_x = [ 0.95 - xD1(1)   ;
                0.95 - xD2(2)   ;
                0.95 - xB2(3)]  ;
        C = [C_x ; Cmax ; Cmin] ;
        Ceq = [] ;
    case('LV')
        Cmin = [-D1 ; -B1 ; -L1 ; -V1 ; -D2 ; -B2 ; -L2 ; -V2] ; %Make sure all flows are positive!
        Cmax = [V1 - Vmax ; V2 - Vmax_2] ;
        C_x = [ 0.95 - xD1(1) ;
                0.95 - xD2(2) ;
                0.95 - xB2(3) ];
        C = [C_x ; Cmax ; Cmin] ;
        Ceq = [] ;
end
        





