slice関数ので作​成した複数のプロット​の同時表示

function func = directivity_soundsource(f,grid_area,interval,number)
%f:周波数 grid_area:代入した値✖代入した値分の領域を表示 interval:表示の細かさ number:２次音源の数

% 定数定義
w = 2*pi*f; %角周波数
lambda = 340/f; %波長
k = 2*pi/lambda;%角波数
K = f/340; %波数
p = 1.293; %空気の密度
c = 340; %音速
A = 1/k^2*p*c;
r = 1; %２次音源を置く円の直径
theta_interval = (2 * pi)/number; %２次音源を置く間隔
xslice = 0;%sliceを入れる座標
yslice = [];
zslice = 0;

%配列定義
main_directivity = [];
grid = [];
directivity = [];
theta_directivity = [];
soundpresure = [];
soundsorece = [];
sum_soundsorce = zeros(1,1,1,1);



% グリッド作成
[X,Y,Z] = meshgrid(-grid_area:interval:grid_area);
soundposition_x = [];
soundposition_y = []; 
soundposition_z = [];
theta = zeros(1,number);

%音源位置計算
for i = 1:number
    theta(1,i) = theta_interval * i;
    for j = 1:3
        if j == 1
            soundposition_x = cat(4,soundposition_x,X - r*cos(theta(1,i)));
        elseif j == 2
            soundposition_y = cat(4,soundposition_y,Y - r*sin(theta(1,i)));
        else
            soundposition_z = cat(4,soundposition_z,Z);
        end
    end
end

%指向性計算
for i = 1:number
    main_directivity = cat(4,main_directivity,[cos(theta(1,i)) sin(theta(1,i)) 0]);
    grid = cat(4,grid,[X(:)-r*cos(theta(1,i)),Y(:)-r*sin(theta(1,i)),Z(:)]);%グリッドの全座標
    grid(:,:,:,i) = normalize(grid(:,:,:,i), 2, 'norm'); % 単位ベクトルに (原点で0で割ってしまうため)
    theta_directivity = cat(4,theta_directivity,reshape(acos(grid(:,:,:,i)*main_directivity(:,:,:,i).'), size(X)));
    directivity = cat(4,directivity,(cos(theta_directivity(:,:,:,i))+1)*0.5);
end

%中心から２次音源までの位相のズレを計算
difference_wave = (r/2)*K;
difference_wave = difference_wave - round(difference_wave); %小数点部分のみ抽出
difference_phase = difference_wave * (2*pi); %位相のズレの計算

%音圧計算
for i = 1:number
    soundsorece = directivity(:,:,:,i).*(1i*w*p*A/(sqrt(soundposition_x(:,:,:,i).^2+soundposition_y(:,:,:,i).^2+soundposition_z(:,:,:,i).^2)).*exp(-1i*k*(sqrt(soundposition_x(:,:,:,i).^2+soundposition_y(:,:,:,i).^2+soundposition_z(:,:,:,i).^2)) - r/2));
%     soundsorece = directivity(:,:,:,i).*exp(1i*k*(sqrt(soundposition_x(:,:,:,i).^2+soundposition_y(:,:,:,i).^2+soundposition_z(:,:,:,i).^2)-r/2));
    soundpresure = cat(4,soundpresure,soundsorece);%グリッドの全座標 
    [value,index] = max(soundpresure(:,:,:,i),[],'all','linear');
    soundpresure(index) = inf;
end

%制御対象の点音源
target_soundsorce = 1i*w*p*A/sqrt(X.^2+Y.^2+Z.^2).*exp(-1i*k*sqrt(X.^2+Y.^2+Z.^2));

for i = 1:number
    sum_soundsorce = soundpresure(:,:,:,i) + sum_soundsorce;
end

%中心の点音源との差分の音圧分布
total_soundpresure1 = abs(sum_soundsorce + target_soundsorce);
%周囲の音源の音圧分布
total_soundpresure2 = abs(sum_soundsorce);
%中心音源のみ
% total_soundpresure = abs(target_soundsorce);

% プロット
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%一つ目にプロットしたい図
slice(X,Y,Z,total_soundpresure1,xslice,yslice,zslice)
title('3次元空間の音圧分布')
xlabel('x[m]') 
ylabel('y[m]') 
zlabel('z[m]')
c = colorbar;
shading flat

%二つ目にプロットしたい図
slice(X,Y,Z,total_soundpresure2,xslice,yslice,zslice)
title('3次元空間の音圧分布')
xlabel('x[m]') 
ylabel('y[m]') 
zlabel('z[m]')
c = colorbar;
shading flat


%以上の二つを同時に画面に出力したいと思っています
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

現在私は特定の音源が作る音圧分布をプロットしようとしています。

音源の配置によってどのような音圧分布になるのかを比べたいと思い、slice関数を使用して複数の状況の同時に画面に表示できないかということを試みたのですが、うまくできず困っています。

slice関数を利用しての複数個のプロットを同時に表示することは可能でしょうか。