diff --git a/doc/manual/source/dynare-misc-commands.rst b/doc/manual/source/dynare-misc-commands.rst
index 7fa7da3eb..b176c3cf3 100644
--- a/doc/manual/source/dynare-misc-commands.rst
+++ b/doc/manual/source/dynare-misc-commands.rst
@@ -53,8 +53,9 @@ Dynare misc commands
 .. command:: generate_trace_plots(CHAIN_NUMBER);
 
     Generates trace plots of the MCMC draws for all estimated
-    parameters and the posterior density in the specified Markov Chain
-    ``CHAIN_NUMBER``.
+    parameters and the posterior density for the specified Markov Chain(s)
+    ``CHAIN_NUMBER``. If ``CHAIN_NUMBER`` is a vector of integers, the trace plots
+    will plot contains separate lines for each chain.
 
     |br|
 
diff --git a/matlab/generate_trace_plots.m b/matlab/generate_trace_plots.m
index a241cc076..409ab22bd 100644
--- a/matlab/generate_trace_plots.m
+++ b/matlab/generate_trace_plots.m
@@ -34,7 +34,7 @@ global M_ options_ estim_params_
 % Get informations about the posterior draws:
 MetropolisFolder = CheckPath('metropolis', M_.dname);
 record=load_last_mh_history_file(MetropolisFolder, M_.fname);
-if chain_number>record.Nblck
+if max(chain_number)>record.Nblck
     error('generate_trace_plots:: chain number is bigger than existing number of chains')
 end
 
diff --git a/matlab/trace_plot.m b/matlab/trace_plot.m
index 8deecb817..369f9adf5 100644
--- a/matlab/trace_plot.m
+++ b/matlab/trace_plot.m
@@ -9,7 +9,8 @@ function trace_plot(options_,M_,estim_params_,type,blck,name1,name2)
 %   type            [string]       'DeepParameter', 'MeasurementError' (for measurement equation error),
 %                                  'StructuralShock' (for structural shock)
 %                                  or 'PosteriorDensity (for posterior density)'
-%   blck            [integer]      Number of the mh chain.
+%   blck            [integer]      Number of the mh chain or chains if a
+%                                  vector
 %   name1           [string]       Object name.
 %   name2           [string]       Object name.
 %
@@ -18,7 +19,7 @@ function trace_plot(options_,M_,estim_params_,type,blck,name1,name2)
 %
 % SPECIAL REQUIREMENTS
 
-% Copyright © 2003-2018 Dynare Team
+% Copyright © 2003-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -62,9 +63,29 @@ TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
 [mh_nblck] = size(record.LastParameters,2);
 clear record;
 
+n_nblocks_to_plot=length(blck);
+
+if n_nblocks_to_plot==1
 % Get all the posterior draws:
 PosteriorDraws = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck);
-
+else
+    PosteriorDraws=NaN(TotalNumberOfMhDraws,n_nblocks_to_plot);
+    save_string='';
+    if options_.TeX
+        title_string_tex='';
+    end
+    for block_iter=1:n_nblocks_to_plot
+        PosteriorDraws(:,block_iter) = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck(block_iter));
+        save_string=[save_string,'_',num2str(blck(block_iter))];
+        if options_.TeX
+            title_string_tex=[title_string_tex, ', ' num2str(blck(block_iter))];
+        end
+    end
+    save_string(1)=[];
+    if options_.TeX
+        title_string_tex(1:2)=[];
+    end
+end
 
 % Plot the posterior draws:
 
@@ -92,30 +113,34 @@ else
     FigureName = ['Trace plot for ' TYPE name1 ' and ' name2];
 end
 
-if options_.mh_nblck>1
-    FigureName = [ FigureName , ' (block number ' int2str(blck)  ').'];
+
+if n_nblocks_to_plot==1
+    if options_.mh_nblck>1
+        FigureName = [ FigureName , ' (block number ' int2str(blck)  ').'];
+    end
+    hh=dyn_figure(options_.nodisplay,'Name',FigureName);
+    plot(1:TotalNumberOfMhDraws,PosteriorDraws,'Color',[.7 .7 .7]);
+
+    % Compute the moving average of the posterior draws:
+    N = options_.trace_plot_ma;
+    MovingAverage = NaN(TotalNumberOfMhDraws,1);
+    first = N+1;
+    last = TotalNumberOfMhDraws-N;
+
+    for t=first:last
+        MovingAverage(t) = mean(PosteriorDraws(t-N:t+N));
+    end
+
+    hold on
+    plot(1:TotalNumberOfMhDraws,MovingAverage,'-k','linewidth',2)
+    hold off
+    axis tight
+    legend({'MCMC draw';[num2str(N) ' period moving average']},'Location','NorthEast')
+else
+    hh=dyn_figure(options_.nodisplay,'Name',FigureName);
+    pp=plot(1:TotalNumberOfMhDraws,PosteriorDraws);
+    legend(pp,strcat(repmat({'Chain '},n_nblocks_to_plot,1),num2str(blck(:))));
 end
-
-hh=dyn_figure(options_.nodisplay,'Name',FigureName);
-plot(1:TotalNumberOfMhDraws,PosteriorDraws,'Color',[.7 .7 .7]);
-
-
-% Compute the moving average of the posterior draws:
-
-N = options_.trace_plot_ma;
-MovingAverage = NaN(TotalNumberOfMhDraws,1);
-first = N+1;
-last = TotalNumberOfMhDraws-N;
-
-for t=first:last
-    MovingAverage(t) = mean(PosteriorDraws(t-N:t+N));
-end
-
-hold on
-plot(1:TotalNumberOfMhDraws,MovingAverage,'-k','linewidth',2)
-hold off
-axis tight
-legend({'MCMC draw';[num2str(N) ' period moving average']},'Location','NorthWest')
 % create subdirectory <dname>/graphs if it doesn't exist
 if ~exist(M_.dname, 'dir')
     mkdir('.',M_.dname);
@@ -130,8 +155,11 @@ if strcmpi(type,'PosteriorDensity')
 else
     plot_name=get_the_name(column,0,M_,estim_params_,options_);
 end
-plot_name=[plot_name,'_blck_',num2str(blck)];
-
+if n_nblocks_to_plot==1
+    plot_name=[plot_name,'_blck_',num2str(blck)];
+else
+    plot_name=[plot_name,'_blck_',save_string];
+end
 dyn_saveas(hh,[M_.dname, filesep, 'graphs', filesep, 'TracePlot_' plot_name],options_.nodisplay,options_.graph_format)
 
 if options_.TeX
@@ -157,10 +185,13 @@ if options_.TeX
             FigureName = ['Trace plot for ' TYPE '$' tex_names{strmatch(name1,base_names,'exact')} '$ and $' tex_names{strmatch(name2,base_names,'exact')} '$'];
         end
     end
-    if options_.mh_nblck>1
-        FigureName = [ FigureName , ' (block number ' int2str(blck)  ').'];
+    if n_nblocks_to_plot==1
+        if options_.mh_nblck>1
+            FigureName = [ FigureName , ' (block number ' int2str(blck)  ').'];
+        end
+    else
+        FigureName = [ FigureName , ' (blocks ' title_string_tex  ').'];
     end
-
     fprintf(fid,'%-s\n','\begin{figure}[H]');
     fprintf(fid,'%-s\n','\centering');
     fprintf(fid,'%-s\n',['  \includegraphics[width=0.8\textwidth]{',[M_.dname, '/graphs/TracePlot_' plot_name],'}\\']);
diff --git a/tests/estimation/slice/fs2000_slice.mod b/tests/estimation/slice/fs2000_slice.mod
index 35a8154e9..ec2f83856 100644
--- a/tests/estimation/slice/fs2000_slice.mod
+++ b/tests/estimation/slice/fs2000_slice.mod
@@ -88,3 +88,7 @@ estimation(order=1,datafile='../fsdat_simul',nobs=192,loglinear,mh_replic=100,mh
 posterior_sampling_method='slice',
 posterior_sampler_options=('rotated',1,'use_mh_covariance_matrix',1)
 );
+
+options_.TeX=1;
+generate_trace_plots(1:2);
+options_.TeX=1;
\ No newline at end of file