function solve(constants)  

%% get input constants
rnat = constants.rnat; 
flag_zlb = constants.flag_zlb;
flag_re = constants.flag_re;
sig_rn = constants.sig_rn;
sig_u = constants.sig_u;
gamma = constants.gamma;

switch flag_re
    case 1
        ygap_em_terminal = NaN;
    case 0


        
        constants_but_re = constants;
        constants_but_re.flag_re = 1;
        try
            load(['workspace/soln_simul_' make_spec_name(constants_but_re) '.mat'],'ygap_em')
            ygap_em_terminal = ygap_em;
            clear ygap_em
        catch err
            error('The simulation of the associated RE solution is not available. This is required for the terminal output gap expectation.')
        end
        clear constants_but_re 
        
end
%% model parameters
% preferences and technology
bbeta       = 1/(1+rnat/400); % discount factor
ies         = 1; % intertemporal elasticity of substitution 
kappa_y     = 0.057; % inflation response to output gap 
kappa_q     =0;
Lambda_pi   = 1; % welfare weight on inflation
Lambda_y    = 0.007; % (relative) welfare weight on output

s_d = 0

CqCy        = 0;
gain        = 0.007;
delta       = 0.03/4;
bbetaU      = 1.0031;

% shocks  
rho_rn      = 0.8;      
rho_xi      = 0.99;
sig_xi      = 0.0165; 

%% globally defined indices for states and controls:
define_idx

%% shock quadrature
switch flag_re
    case 0 
        nn    = [3 3];                          % quadrature nodes for each shock  
        mu    = [0 0];                          % mean vector
        var   = [sig_rn^2 0; 0 sig_xi^2];       % covariance matrix
        [e,w] = qnwnorm(nn,mu,var);             % nodes & weights  
        model.func          = 'model_func';
        
    case 1
        if gamma == 1
        nn    = [3 3];
        mu    = [0 0];
        var   = [sig_rn^2 0; 0 sig_u^2];
        [e,w] = qnwnorm(nn,mu,var);
        model.func          = 'func_AB';
        else
        nn    = [3 3];
        mu    = [0 0];
        var   = [sig_rn^2 0; 0 sig_u^2];
        [e,w] = qnwnorm(nn,mu,var);
        model.func          = 'model_func_re';
        end
end

%% function space 
fspace = set_fspace(constants,bbeta,sig_xi,rho_xi,rho_rn,delta,bbetaU,gamma); 

s_coord = funnode(fspace);
s_nodes = gridmake(s_coord); 
[n_state_space, n_states] = size(s_nodes);
disp(['size state space = ' num2str(n_state_space)])

%% pack model 

% model function
model.params        = {bbeta,delta,ies,kappa_y,kappa_q,Lambda_pi,Lambda_y,s_d,CqCy,ygap_em_terminal,gain,bbetaU,rho_xi,rho_rn,sig_xi,sig_rn,sig_u,flag_zlb,gamma}; % inputs to func in this order
model.discount      = bbeta;
model.e             = e;
model.w             = w;

% model for LQ approximation without ZLB
model_lq = model;
model_lq.params{end} = 0; 

n_controls = 5;

%% solution algorithm parameters
optset('dpsolve','algorithm','newton');             % 'newton' (default) or 'funcit'  
optset('dpsolve','maxit',500);   % 500 (default)
optset('dpsolve','tol',2e-06);                      % 
optset('dpsolve_vmax','maxit',50);                  % 50 (default)
optset('dpsolve_vmax','tol',5e-8);                  % 5e-8 (default)
optset('dpsolve_vmax','maxbacksteps',0);            % 0 (default)
optset('dpsolve_vmax','lcpmethod','minmax');        % 'minmax' (default) or 'smooth'

%% LQ approximation for initial guess
e_init = mu;
x_init = zeros(1,n_controls);
s_init = zeros(1,n_states);  
if flag_re == 0
    logqu_ss = log(1/(1 - bbeta*(1-delta)*1));
    s_init(logqu_idx) = logqu_ss;
    s_init(loglagqu_idx) = logqu_ss; 
end
p_init = zeros(1,n_states);

dpcheck(model,s_init,x_init,e_init);  
tic
[v_lq,x_lq] = lqapprox(model_lq,s_nodes,s_init,x_init,p_init);
fprintf('Time LQ Approximation = %7.2f Seconds\n',toc) 

%% Nonlinear iteration
tic_all = tic;
[coeffs,s_out,v,x_opt] = dpsolve(model,fspace,s_nodes,v_lq,x_lq);    

fprintf('Time Value Function Iteration = %7.2f Seconds\n',toc(tic_all)) 

%% Save results
save(['workspace/soln_' make_spec_name(constants) '.mat'])

