ODE solver - division by zero at time boundaries

Hi

I have a problem when solving my system of ODEs. I simplified it to what I find essential for the problem.

As I see it, the problem is at t = 0, y(1,2) = 0 , whereas dy(2:end,2) = NaN AND at t = p.t_empty, y(1,1) = 0, whereas dy(2:end,1) = NaN

This is however the time interval I am interested in the solution for. It is essentially a system of two tanks. A reaction occurs in tank 1 and the matter is transferred to tank 2 meanwhile.

Do you have any suggestions to how I can solve this?

p.Q = 40;
p.Q_R = 150;
p.V = 20;
p.n = 101;
p.t_empty = p.V/p.Q;

y0 = zeros(p.n*2+2,1);
y0(1:2) = [p.V, 1];

options = odeset('RelTol',1e-12,'AbsTol',1e-12);
[t,y] = ode45('ODE_tank',[0,p.t_empty],y0,options,p);

figure(1)
plot((0:p.n-1),y(end,2:p.n+1))

figure(2)
plot((0:p.n-1),y(end,p.n+3:end))

function dy = ODE_tank(t,y,options,p)

y = reshape(y,[],2);

dy = zeros(size(y));

dy(1,1) = -p.Q;
dy(2,1) = -p.Q_R/y(1,1)*y(2,1);
dy(3:p.n,1) = p.Q_R/y(1,1)*(y(2:p.n-1,1) - y(3:p.n,1));
dy(p.n+1,1) = p.Q_R/y(1,1)*y(p.n,1);

dy(1,2) = p.Q;
dy(2:end,2) = p.Q/y(1,2)*(y(2:end,1)-y(2:end,2));

dy = reshape(dy,[],1);
end

[darova]

• As I see it, the problem is at t = 0, y(1,2) = 0

Can you replace 0 with 1e-3?

[Bastian Andersen]

This works, thank you. How do I accept this as the answer?

Answer 1

[darova]

• As I see it, the problem is at t = 0, y(1,2) = 0

Can you replace 0 with 1e-3?