Using bvp4c to solve a system of n ODEs, where n is an arbitrary number

I am trying to use bvp4c to solve a BVP which consists of a system of n ODEs. Greatly simplified, the ODES are of the form:

y1'=y1*(1-a1)-const*sum(y1,y2,...,yn)

y2'=y2*(1-a2)-const*sum(y1,y2,...,yn)

.

.

yn'=yn*(1-an)-const*sum(y1,y2,...,yn)

where there are n unknown parameters (a1,a2,...,an), so the system of equations is subject to 2n boundary conditions, and const is some constant.

So far I have managed to do this for the simple case of n = 2, using a single m-file with nested functions as per the mat4bvp example in the help file. I have one main function which calls the bvp4c solver, with a nested derivative function which provides the set of ODEs, a nested function which provides the boundary conditions, and an auxiliary function which provides an initial guess for the solution. So for my simple case of n = 2, the derivative function contains two ODEs, the BC function contains 4 BCs, and the initial guess function contains guesses for y1 and y2. This all works great - my derivative function looks something like:

 function dydx = npode(x,y,parameters)
        dydx = [y(1)*(1-parameters(1))-const*(y(1)+y(2))
                y(2)*(1-parameters(2))-const*(y(1)+y(2))];
 end

where parameters is the vector of unknown parameters a1,a2.

However, I am trying to generalise my code to the case of n ODEs, so that the size of the functions called by the solver varies with the number of ODEs to solve. So for example if n = 8, my derivative function dydx would contain 8 entries, my BC function would contain 16 entries, etc. For example, the ith entry of the dydx function would look something like:

 y(i)*(1-parameters(i))-const*(y(1)+y(2)+...+y(i)+....y(n))

Is there an easy way to do this? I have been stuffing around with for loops, int2str and strcat to construct a cell array of dimensions nx1 in which the ith cell contains the string representation of the ith ODE. However, I am unsure where to go from here (can I use eval to do anything useful?) and I am sure there is a much easier and more obvious way. I would be grateful for any suggestions and am happy to provide my m-files if it would help.

Cheers

Greg

Answer 1

[Bjorn Gustavsson]

Wouldnt that just be something like this:

function dydx = npode(x,y,parameters)
      dydx = y.*(1-parameters) - const*sum(y));
end

[Greg]

Thanks Bjorn - for some stupid reason I thought you couldn't use that sort of syntax in a function that was called by the bvp solvers.

[Bjorn Gustavsson]

My pleasure.