Solving coupled odes with ode15s

I have 2 coupled equations

c1*(df/dt)+c2*(df/dz)+c3*(f)+c4*(g)=0

(dg/dt)=c5*f+c6*g

with initial conditions as

f(0,t)=1, g(z,0)=0 and f(z,0)=0

   0<f,g,z,t<1

I converted the first equation into ode by method of lines and want to solve them along with the 2nd equation with ode15s.

This is my code

function brussode(N)
%BRUSSODE  Stiff problem modeling a chemical reaction 

    if nargin < 1
    N = 20;
    end

    tspan = [0; 1];
    y0 = [0; 0];

    options = odeset('Vectorized','on','JPattern',jpattern(N));

    [t,y] = ode15s(@f,tspan,y0,options);

    u = y(:,1:2:end);
    x = (1:N)/(N+1);
    surf(x,t,u);
    view(-40,30);
    xlabel('space');
    ylabel('time');
    zlabel('solution u');
    title(['The Brusselator for N = ' num2str(N)]);
    % --------------------------------------------------------------
    function dydt = f(t,y)
     c1=(0.63*150*(10^-6))/(6.93*(10^-5)*3600);
     c2=1;
     c3=((1-0.63-0.27)*0.0045*150*(10^-6)*5021.8)/(6.93*(10^-5));
     c4=-((1-0.63-0.27)*0.0045*150*(10^-6)*4.62*(10^3))/(1.03*6.93*(10^-5));
     c5=(0.0045*3600*1.03*5021.8)/(4.62*(10^3));
     c6=-0.0045*3600;
     delz=0.05;
     a=c2/(c1*delz)-c3/c1;
     b=-c2/(c1*delz);
     c=c4/(c1*delz);
     dydt = zeros(2*N,size(y,2));      % preallocate dy/dt
     % Evaluate the two components of the function at one edge of 
     % the grid (with edge conditions).
     i = 1;
     dydt(i,:) = 0*a+y(i+1,:).*b-y(i+1,:).*c;
     dydt(i+1,:) = y(i,:).*c5+c6*0;
     % Evaluate the two components of the function at all interior 
     % grid points.
     i = 3:3:2*N-3;
     dydt(i,:) = y(i,:)*a+y(i+1,:)*b-y(i+1,:)*c;
     dydt(i+1,:) = y(i,:).*c5+y(i+1,:).*c6;
     % Evaluate the two components of the function at the other edge 
     % of the grid (with edge conditions).
     i = 2*N-1;
     dydt(i,:) = y(i,:).*a+y(i+1,:).*b-y(i+1,:).*c;
     dydt(i+1,:) = y(i,:).*c5+y(i+1,:).*c6;
     end   % End nested function f
     end   % End function brussode

     function S = jpattern(N)
     B = ones(2*N,5);
     B(2:2:2*N,2) = zeros(N,1);
     B(1:2:2*N-1,4) = zeros(N,1);
     S = spdiags(B,-2:2,2*N,2*N);
     end

But I am getting the error:

Index exceeds matrix dimensions.

Error in brussode/f (line 44)
dydt(i,:) = y(i,:)*a+y(i+1,:)*b-y(i+1,:)*c;

Error in odearguments (line 88)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode15s (line 149)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in brussode (line 13)
[t,y] = ode15s(@f,tspan,y0,options);

How do I correct my code?

Answer 1

[Walter Roberson]

When you use the Vectorized option for odeset, your ode function will be called with an unpredetermined number of columns for y, each corresponding to what would be one y vector if the function had been called without Vectorization turned on. The number of columns could vary between calls to your ode function.

The number of rows for y for your ode function will match length(y0). In your case your y0 is of length 2. But your code in f attempts to access right up to row ((2*N-1)+1) = 2*N . As your y0 is of length 2, your N would have to be 1 for that to work.

The number of rows you output must be the same as the number of rows for the input y. If you have 2*N rows of output you had better make your y0 length 2*N:

y0 = zeros(2*N,1);

Also, in order for your code to be able to work properly, your boundary conditions need to be consistent. You have f(0,t)=1, g(z,0)=0 and f(z,0)=0 so f(0,0)=1 due to the first condition but f(0,0)=0 due to the third condition. Your boundaries are therefore not consistent so after you fix your y0, you will not be able to achieve a useful output.

[Niharika Gupta]

Thanks for your reply Walter. Even if I change the boundary condition to f(0,t)=0, I am getting this graph instead of a parabolic shaped curve. Can you please help me in finding the bug in the code.