Spruce budworm diff.eq,parameter variation using matlab functions

[marellasunny]

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.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]; %This command defines the time interval of interest.
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)
% This command tells MATLAB to use the differential equation solver called ode15s to numerically compute the solution of the differential equation defined by the function f, for time interval tspan, and initial condition u0. The right-hand side of the equation tells MATLAB to store this output in vectors t (for the time points) and u (for the population density).

figure;
drawnow
plot(t,u(:,1)); %This tells MATLAB to plot the solution.
hold on
end
% --------------------------------------------------------------------------

function dudt = f(t,u) %This commands defines the function du/dt = 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)]; %This command inputs the left-hand side of the spruce budworm
%differential equtaion.
end

 

[NemoReally]

I'm not a Matlab user, but I think the technique of adding r and q to the differential equation might work. Essentially, initialize u0 to be a vector with values 0.1 (your current value), r and q. You then define the differential equation f(t,u) to return it's current value as the first element of a 3-vector and 0 for the other 2 elements (r and q being constant throughout each integration). You can then loop through r and q values, something like the following (or however Matlab constructs vectors).

f(t,u) = [u(2).u(1)(1-u(1)/u(3)) + u(3)^2/(1 + u(3)^2) 0 0]

for r = ...
for q = ...
u0 = [0.1 r q]
[t,u] = ode15s(@f,tspan,u0)


The image shows something similar in Mathcad
http://https://www.physicsforums.com/attachment.php?attachmentid=54011&stc=1&d=1355766195
You also might want to check that the equation is valid for your desired r and q ranges ... (eg, what the result if q = 0 ?)