Help correcting my stochastic predator prey code so it will run

I am trying to write a stochastic function for the predator prey model from my differential model, following my textbook however I am unsure if the code below, both the function and the code to run this is correct (I get an error when I run this but not sure how to fix this). Not sure where to tweak the code and if I need to plot three lines or simply two?

Function

function [times,states]=simulate_lotkavolterra(t_final, B, F)

c1 =1; %initail rates

c2 =0.2;

c3 = 0.5;

B = 10; F =2; %initial population of rabbits (B) and foxes (F)

t =0;

states =[B, F];

times =[0];

rates =zeros(3,1) % vector for reaction rates

A = [1 0;

-1 1;

            0  -1];

while t <= t_final;

rates(1) = c1*B;

rates(2) = c2*B*F;

rates(3) = c3*F;

rate = sum(rates); % rate of leaving the state

tau = exprnd(1 / rate); % sojourn time

t = t + tau; % update time

if (rate == 0 t > t_final)

t = t_final;

states = [states; B, F];

times = [times, t];

break;

end

% vector that contains the reaction probabilities:

prob = rates / rate;

index = sample_categorical(prob);

B = B + A(index, 1); % update state

F = F + A(index, 2); % update state

states = [states; B, F];

times = [times, t]; end

Code to run the function

t_final = 60;

%%datapoints = 100; %different step sizes, allows for equally spaced data points

trajectories = 1;

sums = zeros(datapoints, 3);

for i = 1:trajectories

[t_traj, N_traj] = simulate_lotkavolterra(t_final);

[t, N] = time_series(t_traj, N_traj, data points);

sums = sums + N; %sums is zero originally, see line 5 end

averages = sums / trajectories; % scalar multiplication

plot(t, averages(:,1), '-og'); hold on

plot(t, averages(:,2), '-ob');

plot(t, averages(:,3), '-oc');

[Matz Johansson Bergström]

Where is the definition of sample_categorical()?

[Star Strider]

‘(I get an error when I run this but not sure how to fix this)’

Care to elaborate on the error?

[Sarah Macdonald]

So I have managed to fix a few of the problems I was having so my code now reads as below but I don't think that the interaction between predator and prey is correct as I guess it should still be oscillatory. Are my rates (1), (2) and (3) correct?

Function

function [times,states]=simulate_lotkavolterra(t_final)

c1 =1; %initial rates

c2 =0.2;

c3 = 0.5;

B = 10; F =2; %initial population of rabbits (B) and foxes (F)

t =0;

states =[B, F];

times =[0];

rates =zeros(2,1); % vector for reaction rates

A = [1 0; % reaction matrix

-1 1;

            0  -1];

while t <= t_final;

rates(1) = c1*B ;

rates(2) = c2*B*F;

                rates(3) = c3*F;

rate = sum(rates); % rate of leaving the state

tau = exprnd(1 / rate); % sojourn time

t = t + tau; % update time

if (rate == 0 t > t_final)

t = t_final;

states = [states; B, F];

times = [times, t];

break;

end

% vector that contains the reaction probabilities:

prob = rates / rate;

index = sample_categorical(prob);

B = B + A(index, 1); % update state

F = F + A(index, 2); % update state

states = [states; B, F];

times = [times, t]; end

Code to run function

[t,N]=simulate_lotkavolterra (30);

plot (t,N)