function solve(RNAT,flag_zlb,flag_AB,ygapss)
 


rnat = RNAT;
%% model parameters
% preferences and technology
beta       = 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_u      = 0;
sig_u      = 0.154/100; 

sig_rn      = 0.002940; %AB

%% algorithm parameters & function space

spread_psilag = 1;
switch RNAT
    case 1
        tightgrid_lo = -0.28;
        tightgrid_hi = 0.1;
        spread_mulag = 0.5; 
end
 

fspace = fundef(...
            {'spli',4*sig_u/sqrt(1-rho_u^2)*[-1.0 1.0] - (1-rho_u)*kappa_q*s_d,3},... % u
            {'spli',spread_psilag*[-1 1],3},... % psilag
            {'spli',4*sig_rn/sqrt(1-rho_rn^2)*[-1 -0.25:0.05:0.25  1],0},... % rn 
            {'spli',spread_mulag*[-1 1],9}... % mulag
            ); 
        
s_coord = funnode(fspace);
s_nodes = gridmake(s_coord); 
[n_state_space, n_states] = size(s_nodes)

% shock quadrature
nn    = [3 3];                       % quadrature nodes for each shock  
mu    = [0 0];                       % mean vector
var   = [sig_u^2 0; 0 sig_rn^2];     % covariance matrix
[e,w] = qnwnorm(nn,mu,var);          % nodes & weights 


% pack model
model.params        = {beta,ies,kappa_y,kappa_q,delta,Lambda_pi,Lambda_y,s_d,CqCy,ygapss,gain,rho_u,rho_rn,flag_zlb}; % inputs to func in this order
model.discount      = beta;
model.e             = e;
model.w             = w;
    if flag_AB == 1
model.func = 'func_AB';
    else
model.func = 'func';
    end

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

n_controls = 5;

optset('dpsolve','nres',1);                         % add in between nodes for each state; 10 (default) 
optset('dpsolve','algorithm','newton');             % 'newton' (default) or 'funcit' 
optset('dpsolve','maxit',inf);                       % 500 (default)
optset('dpsolve','tol',1e-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 values
e_init = mu;
x_init = zeros(1,n_controls);
s_init = zeros(1,n_states);
s_init(1) = -kappa_q*s_d;
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

[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)

if flag_zlb == 1 
    if flag_AB == 1
        save(['recent_zlb_' num2str(RNAT) '_AB.mat'])
    else
        save(['recent_zlb_' num2str(RNAT) '.mat'])
    end
elseif flag_zlb == 0  
    if flag_AB == 1
        save(['recent_nozlb_' num2str(RNAT) '_AB.mat'])
    else
        save(['recent_nozlb_' num2str(RNAT) '.mat'])
    end
end 