function inspect_residuals(constants)

%% get input constants
rnat = constants.rnat; 
flag_zlb = constants.flag_zlb;
flag_re = constants.flag_re;
 
%% Load results
load(['workspace/soln_' make_spec_name(constants) '.mat'])

max_resid

%% Inspect residuals 
switch flag_re
    case 0 % GUI to go through higher-dimensional array
        xnames={'psilag','mulag','rn','logxi','logqu','laglogqu'}; % leave this unchanged (should correspond to define_idx!)
        perm_order=[3;5;1;2;4;6]; % here, choose how the states should be ordered in the GUI;  
                % the first two are on the axis of a 3D plot
                % the values of the last three are chosen on sliders
        residuals_gui(permute(resid,perm_order),xnames(perm_order),s_coord_resid(perm_order)) 
        
        % collapse to r^n dimension
        figure; plot(s_coord_resid{rn_idx},squeeze(max(abs(resid),[],[psilag_idx mulag_idx logxi_idx logqu_idx loglagqu_idx])),'o-')
        xlabel('$r^n$','Interpreter','Latex'); title('Maximum absolute residual','Interpreter','Latex')
        
    case 1 % no interpolation necessary with just two states  
        xnames={'psilag','mulag','rn'}; % leave this unchanged (should correspond to define_idx!)
        perm_order=[1;3;2]; % here, choose how the states should be ordered in the GUI;  
                % the first two are on the axis of a 3D plot
                % the values of the last three are chosen on sliders
        residuals_gui(permute(resid,perm_order),xnames(perm_order),s_coord_resid(perm_order)) 
        
        % collapse to r^n dimension
        figure; plot(s_coord_resid{rn_idx},squeeze(max(abs(resid),[],[psilag_idx mulag_idx])),'o-')
        xlabel('$r^n$','Interpreter','Latex'); title('Maximum absolute residual','Interpreter','Latex')
        
end  