Несколько графиков в одном рисунке

У меня есть следующий код, и я хочу объединить графики фазового пространства в один рисунок.

Я закодировал функции, но не знаю, как заставить MATLAB объединить их в один рисунок. Как видите, изменяются переменные r, a, b и d. Как мне их объединить?

Я также хотел бы построить векторное поле этих графиков фазового пространства с помощью команды quiver, но она просто не работает.

%function lotkavolterra
% Plots time series and phase space diagrams.
clear all; close all;
t0 = 0;
tf = 20;
N0 = 20;
P0 = 5;

% Original plot
r = 2;
a = 1;
b = 0.2;
d = 1.5;

% Time series plots
lv = @(t,x)(lv_eq(t,x,r,a,b,d));
[t,NP] = ode45(lv,[t0,tf],[N0 P0]);
N = NP(:,1); P = NP(:,2);
figure
plot(t,N,t,P,' --');
axis([0 20 0 50])
xlabel('Time')
ylabel('predator-prey')
title(['r=',num2str(r),', a=',num2str(a),', b=',num2str(b),', d=',num2str(d)]);
saveas(gcf,'predator-prey.png')
legend('prey','predator')

% Phase space plot

figure
quiver(N,P);
axis([0 50 0 10])
%axis tight


% Change variables
r = 2;
a = 1.5;
b = 0.1;
d = 1.5;

%time series plots
lv = @(t,x)(lv_eq(t,x,r,a,b,d));
[t,NP] = ode45(lv,[t0,tf],[N0 P0]);
N = NP(:,1); P = NP(:,2);
figure
plot(t,N,t,P,' --');
axis([0 20 0 50])
xlabel('Time')
ylabel('predator-prey')
title(['r=',num2str(r),', a=',num2str(a),', b=',num2str(b),', d=',num2str(d)]);
saveas(gcf,'predator-prey.png')
legend('prey','predator')


% Phase space plot
figure
plot(N,P);
axis([0 50 0 10])

% Change variables
r = 2;
a = 1;
b = 0.2;
d = 0.5;

% Time series plots
lv = @(t,x)(lv_eq(t,x,r,a,b,d));
[t,NP] = ode45(lv,[t0,tf],[N0 P0]);
N = NP(:,1); P = NP(:,2);
figure
plot(t,N,t,P,' --');
axis([0 20 0 50])
xlabel('Time')
ylabel('predator-prey')
title(['r=',num2str(r),', a=',num2str(a),', b=',num2str(b),', d=',num2str(d)]);
saveas(gcf,'predator-prey.png')
legend('prey','predator')


% Phase space plot
figure
plot(N,P);
axis([0 50 0 10])

% Change variables
r = 0.5;
a = 1;
b = 0.2;
d = 1.5;

% Time series plots
lv = @(t,x)(lv_eq(t,x,r,a,b,d));
[t,NP] = ode45(lv,[t0,tf],[N0 P0]);
N = NP(:,1); P = NP(:,2);
figure
plot(t,N,t,P,' --');
axis([0 20 0 50])
xlabel('Time')
ylabel('predator-prey')
title(['r=',num2str(r),', a=',num2str(a),', b=',num2str(b),', d=',num2str(d)]);
saveas(gcf,'predator-prey.png')
legend('prey','predator')

% Phase space plot
figure
plot(N,P);
axis([0 50 0 10])

% FUNCTION being called from external .m file

%function dx = lv_eq(t,x,r,a,b,d)
%N = x(1);
%P = x(2);
%dN = r*N-a*P*N;
%dP = b*a*P*N-d*P;
%dx =  [dN;dP];