Matlab: adapting myODEfun for the 'Vectorized' propery of odeset

Matlab: adapting myODEfun for the 'Vectorized' propery of odeset

I have a stiff system of coupled ODEs that I am feeding MATLAB's ode15s
solver. It works well, but now I'm trying to optimize the speed of
integration. I am modeling 5 different variables on N different spatial
sites, giving 5N coupled equations. For the moment, N=20 and integration
time is about 25s, but I would like to go to larger values of N.
I used the profiler to see that the vast majority of the time is spent
evaluating myODEfun. I did my best to optimize the code, but that doesn't
change the fact that there is quite a bit going on in the function and
that it is being evaluated ~50,000 times. I read that using the
'Vectorized' property for the ODEfunction can reduce the number of
evaluations needed.
But I don't quite understand what exactly it is that I need to change
about my ODEfun to make it conform to what Matlab wants a 'vectorized'
ODEfun to look like.
From the documentation I see that you can change the example Van der Pol
system from its normal form:
function dydt = vdp1000(t,y)
dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];
to the vectorized form:
function dydt = vdp1000(t,y)
dydt = [y(2,:); 1000*(1-y(1,:).^2).*y(2,:)-y(1,:)];
I don't understand exactly what this new matrix of y is supposed to
represent, and how the size of the second dimension is even defined. I
could almost live with just adding ",:" and not thinking about it, but I
am running into problems because I am already doing some vector operations
in my code.
Here is a simplified example for 2 variables, making 2*N equations. Please
don't try to make sense of the ODEs that are generated here: they don't. I
am talking about the operations that are happening.
function dydt = exampleODEfun(t,y,N)
dydt = zeros(2*N,1);
dTdt = zeros(N,1);
dXdt = zeros(N,1);
T = y(1:N);
X = y(N+1:2*N);
a = [T(2:N).^2 T(2:N) ones(N-1,1)];
b = [3 5 -2];
dTdt(1:N) = 0;
dXdt(1) = 0;
dXdt(2:N) = a*b';
dydt(1:N) = dTdt;
dydt(N+1:2*N) = dXdt;
end
Obviously in the real function a lot more is going on, both for T and X.
As you can see, dXdt(1) is a boundary condition and needs its own
calculations.
Blindly passing odeset 'Vectorized','on' and just adding ",:" to all the
indexes does not work. For example, what size do I need to initialize dTdt
and dXdt to now? What do I do to the ones(N-1,1)? And what do I need to do
to make (a*b') still make sense?
I am using Matlab R2006a.