changing number of variables in ODE

Hi All,

I need to write a function that describes the diffusion of materials in a two-layer sphere (a shell layer and a core layer) into the gas. I use the ode functions in matlab to model the concentrations in gas, shell, and core.

My main function is:

r01 = 8; % initial thickness of the core

r02 = 2; % initial thickness of the shell

d0 = (r01+r02)*2; % initial diameter of the sphere

SA0 = zeros(2,1);

V0 = zeros(2,1);

SA0(1) = pi*d0^2; % initial surface area of the sphere

SA0(2) = pi*(2*r01)^2; % initial surface area the core

V0(1) = pi/6*(d0^3-(2*r01)^3); % initial volume of the shell

V0(2) = pi/6*(2*r01)^3; % initial volume of the core

c0 = [0;V0(1)*1e5;V0(2)*1e5];% initial concentrations

options = odeset('RelTol',1e-10);

[t,c] = ode15s(@myfun,[0:120],c0,options,V0,c0);

subplot(1,2,2)
plot(t,c(:,1),'r.',t,c(:,2),'b.',t,c(:,3),'g.');

The function code is:

function y = myfun(t,c,V0,c0)

V=zeros(2,1);

r=zeros(2,1);

SA=zeros(2,1);

V(1) = V0(1)*c(2)/c0(2);% new volume of the shell

V(2) = V0(2)*c(3)/c0(3);% new volume of the core

d = (6*sum(V)/pi)^(1/3);% new diameter of the sphere

r(2) = (6*V(2)/pi)^(1/3)/2;% new thickness of the core

r(1) = d/2 - r(2);% new thickness of the shell

SA(1) = pi*d^2;% new surface area of the sphere

SA(2) = pi*(2*r(1))^2;% new surface area of the core

y(1) = 2e-6*(1e9-c(1))*SA(1); % gas

y(2) = -2e-6*(1e9-c(1))*SA(1); % shell

y(3) = 0; % core

y = y(:);

subplot(1,2,1)

plot(t, r(1),'ro');hold on;

The thickness of the shell get less than zero at some point, which is not reasonable. I want to force the two layers mix into one layer when the shell thickness reaches below 0.2. Then the materials in the new single layer diffuse into gas at the same rate as the initial shell layer. Any idea how to do that? I understand that this modification will change the number of variables from the initial condition (3 variables into 2 variables). Does ode15s allow this? Thank you so much.

[Walter Roberson]

Please rewrite using @(t,y) myfun(t, y, v0,c0) instead of @myfun,and remove v0 and c0 from the ode15 calling parameters. There are some odd situations that can crop up with the old form that you use.

[Walter Roberson]

You should use a plot function in your options instead of plotting in your myfun