Spruce budworm diff.eq,parameter variation inside the function

I have come across a matlab code for solving the spruce budworm differential equation.But,I would like to solve the same differential equation for a range of parameters(r=0:5,q=0:10).I am having problems trying to define the array of matrices for each loop.How do I define the solution of the differential equation as a array/structure? Please look at the code below:

   % This program solves the differential equation

% du/dt = f(t,u) = ru(1-u/q) + u^2/(1 + u^2)

function spruceco

tspan = [0; 100];

u0 = [.1];

for q=0:10 %I would like to add these here and then pass the values of q,r to the function f%

for r=0:5

[t,u] = ode15s(@f,tspan,u0)

figure;

drawnow

plot(t,u(:,1)); %This tells MATLAB to plot the solution.

hold on

end

% --------------------------------------------------------------------------

function dudt = f(t,u)

q=10;r=2;

%I would like to solve the differential equation for

q=0:10,r=0:5;How?Procedure?

dudt = [r*u(1)*(1-u(1)/q) - u(1)^2/(1+u(1)^2)];

end