function RAPTOR_plot_stationary(timeslice1, timeslice2, varargin)
    set(0, 'defaulttextInterpreter', 'latex');
    set(0, 'defaultLegendInterpreter', 'latex');
    set(0, 'defaultAxesTickLabelInterpreter', 'latex');
    
    if ~isempty(varargin)
      case1 = varargin{1};
      case2 = varargin{2};
    else
      case1 = ' case 1';
      case2 = ' case 2';
    end
    
    figure(1)
    
    hs(1) = subplot(231);
    % plot electron temperature and ec power vs radial coordinate
    
    hs(2) = subplot(232);
    % plot electron density vs radial coordinate
    
    hs(3) = subplot(233);
    % plot chie vs radial coordinate
    
    hs(4) = subplot(234);
    % plot q profile and current drive vs radial coordinate
    
    hs(5) = subplot(235);
    % plot current density (components) vs radial coordinate
    
    hs(6) = subplot(236);
    % plot loop voltage vs radial coordinate
    
    rho1 = timeslice1.rho;
    rho2 = timeslice2.rho;
    
    ax = hs(1);
    plot(ax, rho1, 1e-3*timeslice1.te, 'r--', ...
             rho1, 1e-3*timeslice1.ti, 'r', ...
             rho2, 1e-3*timeslice2.te, 'b--', ...
             rho2, 1e-3*timeslice2.ti, 'b', 'LineWidth', 1.5);
    xlabel(ax, '$\rho$'); 
    legend(ax, {strcat('$T_e$ for', case1), strcat('$T_i$ for', case1), ...
                strcat('$T_e$ for', case2), strcat('$T_i$ for', case2)}, 'Location', 'northeast');
    title(ax, '$T_i$ [$keV$] and $T_e$ [$keV$]');
    ax.FontSize = 14;

    ax = hs(2);
    plot(ax, rho1, 1e-19*timeslice1.ne, 'r--', ...
             rho1, 1e-19*timeslice1.ni, 'r', ...
             rho2, 1e-19*timeslice2.ne, 'b--', ...
             rho2, 1e-19*timeslice2.ni, 'b', 'LineWidth', 1.5);
    xlabel(ax, '$\rho$'); 
    legend(ax, {strcat('$n_e$ for', case1), strcat('$n_i$ for', case1),strcat('$n_e$ for', case2), strcat('$n_i$ for', case2)}, 'Location', 'northeast');
    title(ax, '$n_i$ [$10^{19}m^{-3}$] and $n_e$ [$10^{19}m^{-3}$');
    ax.FontSize = 14;
    
    ax = hs(3);
    plot(ax, rho1, timeslice1.chie, 'r--', ...
             rho1, timeslice1.chii, 'r', ...
             rho1, timeslice1.shear./timeslice1.q, 'r-.', ...
             rho2, timeslice2.chie, 'b--', ... 
             rho2, timeslice2.chii, 'b', ...
             rho2, timeslice2.shear./timeslice2.q, 'b-.', ...
             'LineWidth', 1.5); 
    xlabel(ax, '$\rho$'); 
    legend(ax, {strcat('$\chi_e$ for', case1), strcat('$\chi_i$ for', case1), strcat('$s/q$ for', case1),...
                strcat('$\chi_e$ for', case2), strcat('$\chi_i$ for', case2), strcat('$s/q$ for', case2)},...
                'Location', 'northwest');
    title(ax, '$\chi_i$ [$m^2/s$], $\chi_e$ [$m^2/s$], $s/q$');
    ax.FontSize = 14;
    % rho, timeslice2.dne, 'm', ...
    % rho, 1e2*timeslice2.vne, 'm--', ...
    % strcat('$D_{ne}$ for', case2), strcat('$V_{ne}$ for', case2)}, 
    
    ax = hs(4);
    plot(ax, rho1, timeslice1.q, 'r', ...
             rho1, timeslice1.shear, 'r-.', ...
             rho1, 1e-6*timeslice1.jpar, 'r--', ...
             rho2, timeslice2.q, 'b', ...
             rho2, timeslice2.shear, 'b-.', ...
             rho2, 1e-6*timeslice2.jpar, 'b--', 'LineWidth', 1.5);
    xlabel(ax, '$\rho$'); 
    legend(ax, {strcat('$q$ for', case1), strcat('$s$ for', case1), ...
                strcat('$j_{tot}$ for', case1), ...
                strcat('$q$ for', case2), strcat('$s$ for', case2), ...
                strcat('$j_{tot}$ for', case2)}, 'Location', 'northwest');
    title(ax, '$q$, $s$ and $j_{ec}$ [$MA/m^2$]');
    ax.FontSize = 14;
%     ylim(ax, [0 7]);

    ax = hs(5);
    plot(ax, rho1, 1e-6*timeslice1.pec, 'r-.', ...
             rho1, 1e-6*( timeslice1.pnbe+timeslice1.pnbi ), 'r--', ...
             rho1, 1e-6*( timeslice1.palphae+timeslice1.palphai ), 'r', ...
             rho2, 1e-6*timeslice2.pec, 'b-.', ...
             rho2, 1e-6*( timeslice2.pnbe+timeslice2.pnbi ), 'b--', ...
             rho2, 1e-6*( timeslice2.palphae+timeslice2.palphai ), 'b', 'LineWidth', 1.5);
    xlabel(ax, '$\rho$'); 
    legend(ax, {strcat('$P_{ec}$ for', case1), strcat('$P_{nb}$ for', case1), strcat('$P_{alpha}$ for', case1), ...
                strcat('$P_{ec}$ for', case2), strcat('$P_{nb}$ for', case2), strcat('$P_{alpha}$ for', case2)}, 'Location', 'northwest');
    title(ax, '$p_{ec}$ [$MW/m^2$]');
%     ylim(ax, [0 1]);
    ax.FontSize = 14;
         
    ax = hs(6);
    plot(ax, rho1, 1e-6*timeslice1.joh, 'r', ...
             rho1, 1e-6*timeslice1.jec, 'r--', ...
             rho1, 1e-6*timeslice1.jnb, 'r-.', ...
             rho1, 1e-6*timeslice1.jbs, 'r:', ...
             rho2, 1e-6*timeslice2.joh, 'b', ...
             rho2, 1e-6*timeslice2.jec, 'b--', ...
             rho2, 1e-6*timeslice2.jnb, 'b-.', ...
             rho2, 1e-6*timeslice2.jbs, 'b:', 'LineWidth', 1.5);
    xlabel(ax, '$\rho$');  
    legend(ax, {strcat('$j_{oh}$ for', case1), strcat('$j_{ec}$ for', case1), strcat('$j_{nb}$ for', case1), ...
                strcat('$j_{bs}$ for', case1), strcat('$j_{oh}$ for', case2), ...
                strcat('$j_{ec}$ for', case2), strcat('$j_{nb}$ for', case2), strcat('$j_{bs}$ for', case2)}, 'Location', 'northeast');
    title(ax, 'contributions $j_{tot}$ [$MA/m^2$]');
%     ylim(ax, [-.2 1]);
    ax.FontSize = 14;
    
%     st = suptitle('\textbf{Comparing stationary states}: steady state discharge ($U_{pl}=0V$) and ohmic discharge ($U_{pl}=0.5V$)');
%     st = suptitle('\textbf{Comparing stationary states}: low ($I_{p}=61kA$) and high current solutions ($I_{p}=214kA$)');    
%     st = suptitle('Comparing stationary states);
%     st.FontSize = 16;
    
    set(gcf,'paperpositionmode','auto')

    set(gcf,'position',[100 100 1150 (2/3)*1150])
    
end
