Problems getting answers from ode45

Running this code does not generate any answers, the return matrices only contain NaN values. main file:

clear all;
clc;

%Test Parameters
Kc=0.004;        
Dab=0.006;
Dl=0.0088;

cp0=1;           %Concentration at particle center
c0=ones(2000,1); %Initial conditions 
                 %C(1) to C(1000)=Cp=1, C(1001) to C(2000)=Cf=0

for i=1001:2000; %Initial Cf=0 throughout
    c0(i,1)=0;
end

tspan=linspace(1,200,200);%Time span

[T,C]=ode45(@(t,c) dcdt(t,c,cp0,Kc,Dab,Dl),tspan,c0);

function file:

function dcdt1=dcdt(t,c,cp0,Kc,Dab,Dl)
%Constants
N=2000;         %Total Nodes=2000
r=2.5*10^-3;    %Particle Radius
delr=r/(N/2);   %Node Length for R=1000
z=5;            %Bed Length
delz=z/(N/2);   %Node Length for Z=1000
ez=0.4;         %Bed Voidage
u=2.24*10^-2;   %Fluid Velocity (CO2)

dcdt1=zeros(N,1); %Define system of dcdt

%Node R=0
dcdt1(1)=6*Dab/2*(c(2)-cp0);
%Nodes R=2:999
for i=2:999
    dcdt1(i)=(Dab/delr^2+Dab/(i*delr^2))*c(i+1)...
              -2*(Dab/i*delr^2)*c(i)...
              +(Dab/delr^2-Dab/(i*delr^2))*c(i-1);
end
%Node R=1000
dcdt1(1000)=((-2*Dab/(delr^2))*(4*r^2/delr^2-4*r/delr+1)*c(999)+...
             (2*Dab/(delr^2)*(4*r^2/delr^2-4*r/delr+1)-Kc*(4*r^2/delr^2))*c(1000)+...
              Kc*(4*r^2/delr^2)*c(1001))/(7*r^3/(6*delr^2)-r^2/(2*delr)+r/2-delr/6);
%Nodes Z=1:1000
for i=1001:1999
    dcdt1(i)=(Dl/delz^2-u/2*delz)*c(i+1)...
              -(2/delz^2+(1-ez/ez))*(3*Kc/r)*c(i)+...
              (Dl/delz^2+u/2*delz)*c(i-1)...
              +(1-ez/ez)*(3*Kc/r)*c(1000);
end

There are no syntax errors, but the return variable C is empty. Hope to get any suggestions or answers as to why ode45 is returning this.

[Walter Roberson]

Considering adding an options structure that has function value checking turned on.

[Jan]

The expression "(1-ez/ez)" in the last FOR-loop is at least confusing. Do you mean "(1-ez)/ez"?

In the first FOR-loop you have "Dab/(i*delr^2)" and "Dab/i*delr^2". Is this correct or did you forget the parenthesis in the second expression?

The code can be simplified substantially by using temporary variables e.g. for "delr^2" and "delz^2". Beside the acceleration, this allows for an easier debugging.

for i=1001:2000; c0(i,1)=0; end

Faster and nicer:

c0(1001:2000) = 0;

Answer 1

[Jan]

Examples for a simplification and vectorization - assuming that "(1-ez/ez)" is not a typo:

%Nodes Z=1:1000
d2 = delz ^ 2;
c1 = u / 2 * delz;
dcdt1(1001:1999)=  (Dl/d2 - c1)    * c(1002:2000) ...
                  -6 * Kc/(d2 * r) * c(1001:1999) ...
                  +(Dl/d2 + c1)    * c(1000:1998);

And:

u = 2.24e-2;   % The POWER operator is very expensive!

Note: The leaner and cleaner the code, the easier is finding bugs.

[Red]

Thanks for your reply. Will try to optimize my code and get the logical error out.