Solve System of ODEs with Multiple values of a parameter using vectorization but not looping.

At present I am having a code for plotting solutions of a ode system with multiple initial conditions using vectorization. I want to get the solutions of the same system with a single initial conditon but with multiple values of a parameter (say beta=[0.01;0.02] in code given below). I know how to do it with using for loop but I want to use vectorization instead of looping.

This first code is for multiple initial coditions which run properly.

clc;
clear all;

y0 = 10:200:400
n = length(y0);
p0_all = [50*ones(n,1) y0(:)]';
[t,p] = ode45(@(t,p) lotkasystem(t,p,n),0:.1:15,p0_all);
p = reshape(p,[],n);
nt = length(t);
for k = 1:n
    plot(p(1:nt,k),p(nt+1:2*nt,k))
    hold on
end
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
hold off


function dpdt = lotkasystem(t,p,n)
%LOTKA  Lotka-Volterra predator-prey model for system of inputs p.
delta = 0.02;
beta = 0.01;

% Change the size of p to be: Number of equations-by-number of initial
% conditions.
p = reshape(p,[],n);

% Write equations in vectorized form.
dpdt = [p(1,:) .* (1 - beta*p(2,:)); 
        p(2,:) .* (-1 + delta*p(1,:))];

% Linearize output.
dpdt = dpdt(:);
end

This second code is for single initial condition and multiple values of beta which is giving error. Please help to rslove this.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;
clear all;

y0 = 10
beta=[0.01; 0.02]';
n = length(beta);

[t,p] = ode45(@(t,p) lotkasystem(t,p,n),0:.1:15,[10 10]);
p = reshape(p,[],n);
nt = length(t);
for k = 1:n
    plot(p(1:nt,k),p(nt+1:2*nt,k))
    hold on
end
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
hold off


function dpdt = lotkasystem(t,p,n)
%LOTKA  Lotka-Volterra predator-prey model for system of inputs p.
delta = 0.02;
beta=[0.0; 0.02]'

% Change the size of p to be: Number of equations-by-number of initial
% conditions.
p = reshape(p,[],n);

% Write equations in vectorized form.
dpdt = [p(1,:) .* (1 - beta*p(2,:)); 
        p(2,:) .* (-1 + delta*p(1,:))];

% Linearize output.
dpdt = dpdt(:);
end

Answer 1

[VBBV]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;
clear all;

y0 = 10:200:400

y0 = 1×2

10 210

n = length(y0);
p0_all = [50*ones(n,1) y0(:)]';
beta=[0.0; 0.02]';
n = length(beta);

[t,p] = ode45(@(t,p) lotkasystem(t,p,n),0:.1:15,p0_all);
p = reshape(p,[],n);
nt = length(t);
for k = 1:n
    plot(p(1:nt,k),p(nt+1:2*nt,k))
    hold on
end
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
hold off

function dpdt = lotkasystem(t,p,n)
%LOTKA  Lotka-Volterra predator-prey model for system of inputs p.
delta = 0.02;
beta=[0.01; 0.02]';  % change 0.0 with small value as defined before. 

% Change the size of p to be: Number of equations-by-number of initial
% conditions.
p = reshape(p,[],n);

% Write equations in vectorized form.
dpdt = [p(1,:) .* (1 - beta.*p(2,:));   % 
        p(2,:) .* (-1 + delta*p(1,:))];

% Linearize output.
dpdt = dpdt(:);
end

[VBBV]

dpdt returns a 4x1 vector. Give beta a small value like in earlier case instead of 0 which result in NaN

[VBBV]

clc;
clear all;

y0 = 10:200:400

y0 = 1×2

10 210

n = length(y0);
p0_all = [50*ones(n,1) y0(:)]';
[t,p] = ode45(@(t,p) lotkasystem(t,p,n),0:.1:15,p0_all);
p = reshape(p,[],n);
nt = length(t);
for k = 1:n
    plot(p(1:nt,k),p(nt+1:2*nt,k))
    hold on
end
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
hold off

function dpdt = lotkasystem(t,p,n)
%LOTKA  Lotka-Volterra predator-prey model for system of inputs p.
delta = 0.02;
beta = 0.01;

% Change the size of p to be: Number of equations-by-number of initial
% conditions.
p = reshape(p,[],n);

% Write equations in vectorized form.
dpdt = [p(1,:) .* (1 - beta*p(2,:)); 
        p(2,:) .* (-1 + delta*p(1,:))];

% Linearize output.
dpdt = dpdt(:);
end

[VBBV]

this is with single value of beta and you can see the difference in intial conditions too

[SAJAN Phutela]

I edited the question by taking non-zero beta. Thanks a lot for your help.