Error from ode15s

Received this error: Error in ode15s (line 153)

odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Below is my code:

clc

psi = 0.63;

psip = 0.27;

L = 150e-6;

Ts = 60;

k = 0.00158;

q0 = 4.62e3;

C0 =1.03;

gama = 5021.8;

Jw = 8.5e-6;

a = psi*L/Jw/Ts;

b = (1-psi-psip)*k*L/Jw;

e = q0./C0;

d = k*Ts*C0./q0;

%TIME SPAN

tt = 0:100:2400;

t = tt/Ts;

%STEP LENGTH

Nz = 51;

zz= linspace(0, L, Nz);

z = zz./L;

dz = z(2) -z(1);

%IC

IC_A = zeros(1, Nz);

IC_B = zeros(1, Nz);

IC = [IC_A, IC_B];

%Solve the pde using ode15s

[t, c] = ode15s(@f, t, IC, [], Nz, C0, a, b, e, d, dz, gama);

% Extract solutions

C = c(:, 1:Nz);

g = c(:, Nz+1 : 2*Nz);

% Reinput BC

C(:, 1) = 1;

C(:, end) = (4*C(:, Nz-1)-C(:, Nz-2))./3;

%plot

figure(1)

plot(t, C(:, end), '--')

ylim([0 1.2])

function dcqdt = f(t, c, Nz, C0, a, b, e, d, dz, gama)

dcqdt = zeros(length(c), 1);

dcdt = zeros(Nz, 1);

dqdt = zeros(Nz, 1);

%define value

C = c(1:Nz);

g = c(Nz+1 : 2*Nz);

%BC

C(1) =1;

C(end) = (4*C(Nz-1)-C(Nz-2))./3;

%interior for-loop

for i=2:Nz-1

dcdz(i) = C(i+1) - C(i-1)./2./dz; % Using FDA

dcdt(i) = (-1./a).*(dcdz(i) + b.*(gama.*C(i)-q(i).*e));

dqdt(i) = d.*(gama.*C(i)-q(i).*e);

end

dcqdt = [dcdt; dqdt];

end

Answer 1

[Torsten]

psi = 0.63;
psip = 0.27;
L = 150e-6;
Ts = 60;
k = 0.00158;
q0 = 4.62e3;
C0 =1.03; 
gama = 5021.8;
Jw = 8.5e-6;
a = psi*L/Jw/Ts;
b = (1-psi-psip)*k*L/Jw;
e = q0./C0;
d = k*Ts*C0./q0;
%TIME SPAN
tt = 0:100:24000;
t = tt/Ts;
%STEP LENGTH
Nz = 51;
zz= linspace(0, L, Nz);
z = zz./L;
dz = z(2) -z(1);
%IC
IC_A = zeros(1, Nz);
IC_B = zeros(1, Nz);
IC = [IC_A, IC_B];
%Solve the pde using ode15s
[t, c] = ode15s(@(t,y)f(t, y, Nz, C0, a, b, e, d, dz, gama), t, IC);
% Extract solutions
C = c(:, 1:Nz);
g = c(:, Nz+1 : 2*Nz); 
% Reinput BC
C(:, 1) = 1;
C(:, end) = (4*C(:, Nz-1)-C(:, Nz-2))./3;
%plot
figure(1)
plot(t, C(:, end), '--')

%ylim([0 1.2])
function dcqdt = f(t, c, Nz, C0, a, b, e, d, dz, gama)
dcqdt = zeros(length(c), 1);
dcdt = zeros(Nz, 1);
dqdt = zeros(Nz, 1);
%define value
C = c(1:Nz);
q = c(Nz+1 : 2*Nz); 
%BC
C(1) =1;
C(end) = (4*C(Nz-1)-C(Nz-2))./3;
%interior for-loop 
for i=2:Nz-1
    dcdz(i) = (C(i+1) - C(i-1))./2./dz; % Using FDA
    dcdt(i) = (-1./a).*(dcdz(i) + b.*(gama.*C(i)-q(i).*e));
    dqdt(i) = d.*(gama.*C(i)-q(i).*e);
end
dcqdt = [dcdt; dqdt];
end

[University Glasgow]

Thank very much.

[University Glasgow]

Hi Torsten,

Why is the plot different from what I'm expecting? I have attached the expected plot.

[Torsten]

Why is the plot different from what I'm expecting?

Set

tt = 0:100:24000

instead of

tt = 0:100:2400

(see above)

[University Glasgow]

Okay, thank you