2D PDEs in Cylindrical Coordinates, odes15s, indices compatibility issue

Help will be highly appriciated.

In code line 60 - 68. If I choose mf=1, code is working perectly. But if I put mf=2 it shows some errors which are connected to pde_2.

% ODE integration

reltol=1.0e-04; abstol=1.0e-04;

options=odeset('RelTol',reltol,'AbsTol',abstol);

mf=1;%% If we put mf=2 Code is not working 11/01/2024 at line 62.

if(mf==1)

[t,y]=ode15s(@pde_1,tout,y0,options);

end

if(mf==2)

[t,y]=ode15s(@pde_2,tout,y0,options);

end

%
% Clear previous files
clear all
clc
%
% Global area
global nr nz dr dz drs dzs r z Dc Dt ca Tk cae Tke h k E R rk0 v rho Cp Tw dH ncall
%
% Model parameters
ca0=0.0;
cae=0.01;
Tk0=305.0;
Tke=305.0;
Tw=355.0;
r0=2.0;
zl=100.0;
v=1.0;
Dc=0.1;
Dt=0.1;
k=0.01;
h=0.01;
rho=1.0;
Cp=0.5;
rk0=1.5e+09;
dH=-10000.0;
E=15000.0;
R=1.987;
% Grid in axial direction
nz=20;
dz=zl/nz;
for i=1:nz
z(i)=i*dz;
end
%
% Grid in radial direction
nr=7;
dr=r0/(nr-1);
for j=1:nr
r(j)=(j-1)*dr;
end
drs=dr^2;
%
% Independent variable for ODE integration
tf=200.0;
tout=[0.0:50.0:tf]';
nout=5;
ncall=0;
%
% Initial condition
for i=1:nz
for j=1:nr
ca(i,j)=ca0;
Tk(i,j)=Tk0;
y0((i-1)*nr+j)=ca(i,j);
y0((i-1)*nr+j+nz*nr)=Tk(i,j);
end
end
%
% ODE integration
reltol=1.0e-04; abstol=1.0e-04;
options=odeset('RelTol',reltol,'AbsTol',abstol);
mf=2;%% If we put mf=2 Code is not working 20/9/2023
if(mf==1)
    [t,y]=ode15s(@pde_1,tout,y0,options);
end
if(mf==2)
    [t,y]=ode15s(@pde_2,tout,y0,options);
end

r4fdx =

     []

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Error in solution>dss004 (line 672)
ux( 1)=r4fdx*...

Error in solution>pde_2 (line 206)
car1d=dss004(0.0,r0,nr,ca1d)

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

Error in ode15s (line 148)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

%
% 1D to 2D matrices
for it=1:nout
for i=1:nz
for j=1:nr
ca(it,i,j)=y(it,(i-1)*nr+j);
Tk(it,i,j)=y(it,(i-1)*nr+j+nz*nr);
end
end
end
%
% Display a heading and centerline output
fprintf('\n nr = %2d nz = %2d\n',nr,nz);
for it=1:nout
fprintf('\n t = %4.1f\n',t(it));
for i=1:nz
fprintf(' z = %5.1f ca(r=0,z,t) = %8.5f  Tk(r=0,z,t) = %8.2f\n',z(i),ca(it,i,1),Tk(it,i,1));
end
end
fprintf('\n ncall = %5d\n',ncall);
%
% Parametric plots
%
% Axial profiles
subplot(2,2,1)
plot(z,ca(2,:,1)); axis([0 zl 0 0.01]);
title('ca(r=0,z,t=50)'); xlabel('z');
ylabel('ca(r=0,z,t=50)')
subplot(2,2,2)
plot(z,Tk(2,:,1)); axis([0 zl 305 450]);
title('Tk(r=0,z,t=50)'); xlabel('z');
ylabel('Tk(r=0,z,t=50)')
subplot(2,2,3)
plot(z,ca(3,:,1)); axis([0 zl 0 0.01]);
title('ca(r=0,z,t=100)'); xlabel('z');
ylabel('ca(r=0,z,t=100)')
subplot(2,2,4)
plot(z,Tk(3,:,1)); axis([0 zl 305 450]);
title('Tk(r=0,z,t=100)'); xlabel('z');
ylabel('Tk(r=0,z,t=100)')
figure(2);
subplot(2,2,1)
plot(z,ca(4,:,1)); axis([0 zl 0 0.01]);
title('ca(r=0,z,t=150)'); xlabel('z');
ylabel('ca(r=0,z,t=150)')
subplot(2,2,2)
plot(z,Tk(4,:,1)); axis([0 zl 305 450]);
title('Tk(r=0,z,t=150)'); xlabel('z');
ylabel('Tk(r=0,z,t=150)')
subplot(2,2,3)
plot(z,ca(5,:,1)); axis([0 zl 0 0.01]);
title('ca(r=0,z,t=200)'); xlabel('z');
ylabel('ca(r=0,z,t=200)')
subplot(2,2,4)
plot(z,Tk(5,:,1)); axis([0 zl 305 450]);
title('Tk(r=0,z,t=200)'); xlabel('z');
ylabel('Tk(r=0,z,t=200)')
% Radial profiles (at t = tf)
for i=1:nz,
for j=1:nr
ca_rad(i,j)=ca(5,i,j);
Tk_rad(i,j)=Tk(5,i,j);
end
end
figure(3);
subplot(2,2,1)
plot(r,ca_rad(1,:)); axis([0 r0 0.0 0.01]);
title('ca(r,z=5,t=200)'); xlabel('r');
ylabel('ca(r,z=5,t=200)')
subplot(2,2,2)
plot(r,ca_rad(7,:)); axis([0 r0 0.0 0.01]);
title('ca(r,z=35,t=200)'); xlabel('r');
ylabel('ca(r,z=35,t=200)')
subplot(2,2,3)
plot(r,ca_rad(14,:)); axis([0 r0 0.0 0.01]);
title('ca(r,z=70,t=200)'); xlabel('r');
ylabel('ca(r,z=70,t=200)')
subplot(2,2,4)
plot(r,ca_rad(20,:)); axis([0 r0 0.0 0.01]);
%
title('ca(r,z=100,t=200)'); xlabel('r');
ylabel('ca(r,z=100,t=200)')
figure(4);
subplot(2,2,1)
plot(r,Tk_rad(1,:)); axis([0 r0 305 450]);
title('Tk(r,z=5,t=200)'); xlabel('r');
ylabel('Tk(r,z=5,t=200)')
subplot(2,2,2)
plot(r,Tk_rad(7,:)); axis([0 r0 305 450]);
title('Tk(r,z=35,t=200)'); xlabel('r');
ylabel('Tk(r,z=35,t=200)')
subplot(2,2,3)
plot(r,Tk_rad(14,:)); axis([0 r0 305 450]);
title('Tk(r,z=70,t=200)'); xlabel('r');
ylabel('Tk(r,z=70,t=200)')
subplot(2,2,4)
plot(r,Tk_rad(20,:)); axis([0 r0 305 450]);
title('Tk(r,z=100,t=200)'); xlabel('r');
ylabel('Tk(r,z=100,t=200)')
% 3D plotting
figure()
surf(r,z,Tk_rad)
axis([0 r0 0 zl 0 500]);
xlabel('r - radial distance');
ylabel('z - axial distance');
zlabel('Tk - reactant temperature');
title([{'Surface & contour plot of reactant temperature,Tk(r,z),'},{ 'at time t=200.0'}]);
view(-130,62);
figure()
surf(r,z,ca_rad)
axis([0 r0 0 zl 0 0.01]);
xlabel('r - radial distance');
ylabel('z - axial distance');
zlabel('ca - reactant concentration');
title([{'Surface & contour plot of reactant concentration,ca(r,z),'},{ 'at time t=200.0'}]);
view(-130,62);
% rotate3d on
function yt=pde_2(t,y)
%
% Global area
global nr nz dr dz drs dzs r r0 z zl Dc Dt ca Tk cae Tke h k E R rk0 v rho Cp Tw dH ncall
%
% 1D to 2D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
ca(i,j)=y(ij);
Tk(i,j)=y(ij+nr*nz);
end
end
%
% Step through the grid points in r and z
for i=1:nz
ca1d=ca(i,:);
Tk1d=Tk(i,:);
%
% car, Tkr
car1d=dss004(0.0,r0,nr,ca1d)
car(i,:)=car1d;
car(i,1)= 0.0;
car(i,nr)=0.0;
Tkr1d=dss004(0.0,r0,nr,Tk1d)
Tkr(i,:)=Tkr1d;
Tkr(i,1)= 0.0;
Tkr(i,nr)=(h/k)*(Tw-Tk(i,nr))
%
% carr, Tkrr
car1d( 1)=0.0;
car1d(nr)=0.0;
nl=2;
nu=2;
carr1d=dss044(0.0,r0,nr,ca1d,car1d,nl,nu)
carr(i,:)=carr1d;
Tkr1d( 1)=0.0;
Tkr1d(nr)=(h/k)*(Tw-Tk1d(nr));
nl=2;
nu=2;
Tkrr1d=dss044(0.0,r0,nr,Tk1d,Tkr1d,nl,nu)
Tkrr(i,:)=Tkrr1d;
for j=1:nr
%
% (1/r)*car, (1/r)*Tkr
if(j~=1)
car(i,j)=(1.0/r(j))*car(i,j);
Tkr(i,j)=(1.0/r(j))*Tkr(i,j);
end
if(j==1)
    carr(i,j)=2.0*carr(i,j);
end
if(j==1)
    Tkrr(i,j)=2.0*Tkrr(i,j);
end
%
% caz, Tkz
if(i==1)
caz(i,j)=(ca(i,j)-cae)/dz;
Tkz(i,j)=(Tk(i,j)-Tke)/dz;
else
caz(i,j)=(ca(i,j)-ca(i-1,j))/dz;
Tkz(i,j)=(Tk(i,j)-Tk(i-1,j))/dz;
end
%
% PDEs
rk=rk0*exp(-E/(R*Tk(i,j)))*ca(i,j)^2;
cat(i,j)=Dc*(carr(i,j)+car(i,j))-v*caz(i,j)-rk;
Tkt(i,j)=Dt*(Tkrr(i,j)+Tkr(i,j))-v*Tkz(i,j)-dH/(rho*Cp)*rk;
end
end
%
% 2D to 1D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
yt(ij)=cat(i,j);
yt(ij+nr*nz)=Tkt(i,j);
end
end
%
% Transpose and count
yt=yt';
ncall=ncall+1;
end
function yt=pde_1(t,y)
%
% Global area
global nr nz dr dz drs dzs r z Dc Dt ca Tk cae Tke h k E R rk0 v rho Cp Tw dH ncall
%
% 1D to 2D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
ca(i,j)=y(ij);
Tk(i,j)=y(ij+nr*nz);
end
end
%
% Step through the grid points in r and z
for i=1:nz
for j=1:nr
%
% (1/r)*car, (1/r)*Tkr
if(j==1)
car(i,j)=2.0*(ca(i,j+1)-ca(i,j))/drs;
Tkr(i,j)=2.0*(Tk(i,j+1)-Tk(i,j))/drs;
elseif(j==nr)
car(i,j)=0.0;
Tkr(i,j)=(1.0/r(j))*(h/k)*(Tw-Tk(i,j));
else
car(i,j)=(1.0/r(j))*(ca(i,j+1)-ca(i,j-1))/(2.0*dr);
Tkr(i,j)=(1.0/r(j))*(Tk(i,j+1)-Tk(i,j-1))/(2.0*dr);
end
%
% carr, Tkrr
if(j==1)
carr(i,j)=2.0*(ca(i,j+1)-ca(i,j))/drs;
Tkrr(i,j)=2.0*(Tk(i,j+1)-Tk(i,j))/drs;
elseif(j==nr)
carr(i,j)=2.0*(ca(i,j-1)-ca(i,j))/drs;
Tkf=Tk(i,j-1)+2.0*dr*h/k*(Tw-Tk(i,j));
Tkrr(i,j)=(Tkf-2.0*Tk(i,j)+Tk(i,j-1))/drs;
else
carr(i,j)=(ca(i,j+1)-2.0*ca(i,j)+ca(i,j-1))/drs;
Tkrr(i,j)=(Tk(i,j+1)-2.0*Tk(i,j)+Tk(i,j-1))/drs;
end
%
% caz, Tkz
if(i==1)
caz(i,j)=(ca(i,j)-cae)/dz;
Tkz(i,j)=(Tk(i,j)-Tke)/dz;
else
caz(i,j)=(ca(i,j)-ca(i-1,j))/dz;
Tkz(i,j)=(Tk(i,j)-Tk(i-1,j))/dz;
end
%
% PDEs
rk=rk0*exp(-E/(R*Tk(i,j)))*ca(i,j)^2;
cat(i,j)=Dc*(carr(i,j)+car(i,j))-v*caz(i,j)-rk;
Tkt(i,j)=Dt*(Tkrr(i,j)+Tkr(i,j))-v*Tkz(i,j)-dH/(rho*Cp)*rk;
end
end
%
% 2D to 1D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
yt(ij)=cat(i,j);
yt(ij+nr*nz)=Tkt(i,j);
end
end
%
% Transpose and count
yt=yt';
ncall=ncall+1;
end
function [ux]=dss004(xl,xu,n,u)

% From Chapter 9, "Fitzhugh-Nagumo Equation", of the book: 
% Graham W Griffiths and William E Schiesser (2011),
% Traveling Wave Analysis of Partial Differential Equations
% - Numerical and Analytical Methods with Matlab and Maple,
% Academic Press (ISBN: 9780123846525).   
%
%  Function dss004 computes the first derivative, u , of a
%                                                  x
%  variable u over the spatial domain xl le x le xu from classical
%  five-point, fourth-order finite difference approximations
%
%  Argument list
%
%     xl      Lower boundary value of x (input)
%
%     xu      Upper boundary value of x (input)
%
%     n       Number of grid points in the x domain including the
%             boundary points (input)
%
%     u       One-dimensional array containing the values of u at
%             the n grid point points for which the derivative is
%             to be computed (input)
%
%     ux      One-dimensional array containing the numerical
%             values of the derivatives of u at the n grid points
%             (output)
%
%  The mathematical details of the following Taylor series (or
%  polynomials) are given in routine dss002.
%
%  Five-point formulas
%
%     (1)  Left end, point i = 1
%
%                                   2            3            4
%  a(u2 = u1 + u1  ( dx) + u1  ( dx)  + u1  ( dx)  + u1  ( dx)
%                x   1f      2x  2f       3x  3f       4x  4f
%
%                      5            6            7
%           + u1  ( dx)  + u1  ( dx)  + u1  ( dx)  + ...)
%               5x  5f       6x  6f       7x  7f
%
%                                   2            3            4
%  b(u3 = u1 + u1  (2dx) + u1  (2dx)  + u1  (2dx)  + u1  (2dx)
%                x   1f      2x  2f       3x  3f       4x  4f
%
%                      5            6            7
%           + u1  (2dx)  + u1  (2dx)  + u1  (2dx)  + ...)
%               5x  5f       6x  6f       7x  7f
%
%                                   2            3            4
%  c(u4 = u1 + u1  (3dx) + u1  (3dx)  + u1  (3dx)  + u1  (3dx)
%                x   1f      2x  2f       3x  3f       4x  4f
%
%                      5            6            7
%           + u1  (3dx)  + u1  (3dx)  + u1  (3dx)  + ...)
%               5x  5f       6x  6f       7x  7f
%
%                                   2            3            4
%  d(u5 = u1 + u1  (4dx) + u1  (4dx)  + u1  (4dx)  + u1  (4dx)
%                x   1f      2x  2f       3x  3f       4x  4f
%
%                      5            6            7
%           + u1  (4dx)  + u1  (4dx)  + u1  (4dx)  + ...)
%               5x  5f       6x  6f       7x  7f
%
%  Constants a, b, c and d are selected so that the coefficients
%  of the u1  terms sum to one and the coefficients of the u1  ,
%           x                                                2x
%  u1   and u1   terms sum to zero
%    3x       4x
%
%  a +   2b +   3c +   4d = 1
%
%  a +   4b +   9c +  16d = 0
%
%  a +   8b +  27c +  64d = 0
%
%  a +  16b +  81c + 256d = 0
%
%  Simultaneous solution for a, b, c and d followed by the solu-
%  tion of the preceding Taylor series, truncated after the u
%                                                            4x
%  terms, for u1  gives the following five-point approximation
%               x
%                                                         4
%  u1  = (1/12dx)(-25u1 + 48u2 - 36u3 + 16u4 - 3u5) + O(dx )   (1)
%    x
%
%     (2)  Interior point, i = 2
%
%                                   2            3            4
%  a(u1 = u2 + u2  (-dx) + u2  (-dx)  + u2  (-dx)  + u2  (-dx)
%                x   1f      2x  2f       3x  3f       4x  4f
%
%                      5            6            7
%           + u2  (-dx)  + u2  (-dx)  + u2  (-dx)  + ...)
%               5x  5f       6x  6f       7x  7f
%
%                                   2            3            4
%  b(u3 = u2 + u2  ( dx) + u2  ( dx)  + u2  ( dx)  + u2  ( dx)
%                x   1f      2x  2f       3x  3f       4x  4f
%
%                      5            6            7
%           + u2  ( dx)  + u2  ( dx)  + u2  ( dx)  + ...)
%               5x  5f       6x  6f       7x  7f
%
%                                   2            3            4
%  c(u4 = u2 + u2  (2dx) + u2  (2dx)  + u2  (2dx)  + u2  (2dx)
%                x   1f      2x  2f       3x  3f       4x  4f
%
%                      5            6            7
%           + u2  (2dx)  + u2  (2dx)  + u2  (2dx)  + ...)
%               5x  5f       6x  6f       7x  7f
%
%                                   2            3            4
%  d(u5 = u2 + u2  (3dx) + u2  (3dx)  + u2  (3dx)  + u2  (3dx)
%                x   1f      2x  2f       3x  3f       4x  4f
%
%                      5            6            7
%           + u2  (3dx)  + u2  (3dx)  + u2  (3dx)  + ...)
%               5x  5f       6x  6f       7x  7f
%
%  -a +   b +  2c +  3d = 1
%
%   a +   b +  4c +  9d = 0
%
%  -a +   b +  8c + 27d = 0
%
%   a +   b + 16c + 81d = 0
%
%  Simultaneous solution for a, b, c and d followed by the solu-
%  tion of the preceding Taylor series, truncated after the u
%                                                            4x
%  terms, for u1  gives the following five-point approximation
%               x
%                                                        4
%  u2  = (1/12dx)(-3u1 - 10u2 + 18u3 -  6u4 +  u5) + O(dx )    (2)
%    x
%
%     (3)  Interior point i, i ne 2, n-1
%
%                                        2             3
%  a(ui-2 = ui + ui  (-2dx)  + ui  (-2dx)  + ui  (-2dx)
%                  x    1f       2x   2f       3x   3f
%
%                          4             5             6
%              + ui  (-2dx)  + ui  (-2dx)  + ui  (-2dx)  + ...)
%                  4x   4f       5x   5f       6x   6f
%
%                                        2             3
%  b(ui-1 = ui + ui  ( -dx)  + ui  ( -dx)  + ui  ( -dx)
%                  x    1f       2x   2f       3x   3f
%
%                          4             5             6
%              + ui  ( -dx)  + ui  ( -dx)  + ui  ( -dx)  + ...)
%                  4x   4f       5x   5f       6x   6f
%
%                                        2             3
%  c(ui+1 = ui + ui  (  dx)  + ui  (  dx)  + ui  (  dx)
%                  x    1f       2x   2f       3x   3f
%
%                          4             5             6
%              + ui  (  dx)  + ui  (  dx)  + ui  (  dx)  + ...)
%                  4x   4f       5x   5f       6x   6f
%
%                                        2             3
%  d(ui+2 = ui + ui  ( 2dx)  + ui  ( 2dx)  + ui  ( 2dx)
%                  x    1f       2x   2f       3x   3f
%
%                          4             5             6
%              + ui  ( 2dx)  + ui  ( 2dx)  + ui  ( 2dx)  + ...)
%                  4x   4f       5x   5f       6x   6f
%
%   -2a -   b +   c +  2d = 1
%
%    4a +   b +   c +  4d = 0
%
%   -8a -   b +   c +  8d = 0
%
%   16a +   b +   c + 16d = 0
%
%  Simultaneous solution for a, b, c and d followed by the solu-
%  tion of the preceding Taylor series, truncated after the u
%                                                            4x
%  terms, for u1  gives the following five-point approximation
%               x
%                                                          4
%  ui  = (1/12dx)(ui-2 - 8ui-1 + 0ui + 8ui+1 - ui+2) + O(dx )  (3)
%    x
%
%     (4)  Interior point, i = n-1
%
%                                              2               3
%  a(un-4 = un-1 + un-1  (-3dx)  + un-1  (-3dx)  + un-1  (-3dx)
%                      x    1f         2x   2f         3x   3f
%
%                       4               5               6
%         + un-1  (-3dx)  + un-1  (-3dx)  + un-1  (-3dx)  + ...
%               4x   4f         5x   5f         6x   6f
%
%                                              2               3
%  b(un-3 = un-1 + un-1  (-2dx)  + un-1  (-2dx)  + un-1  (-2dx)
%                      x    1f         2x   2f         3x   3f
%
%                       4               5               6
%         + un-1  (-2dx)  + un-1  (-2dx)  + un-1  (-2dx)  + ...
%               4x   4f         5x   5f         6x   6f
%
%                                              2               3
%  c(un-2 = un-1 + un-1  ( -dx)  + un-1  (- -x)  + un-1  ( -dx)
%                      x    1f         2x   2f         3x   3f
%
%                       4               5               6
%         + un-1  ( -dx)  + un-1  ( -dx)  + un-1  ( -dx)  + ...
%               4x   4f         5x   5f         6x   6f
%
%                                              2               3
%  d(un   = un-1 + un-1  (  dx)  + un-1  (  dx)  + un-1  (  dx)
%                      x    1f         2x   2f         3x   3f
%
%                       4               5               6
%         + un-1  (  dx)  + un-1  (  dx)  + un-1  (  dx)  + ...
%               4x   4f         5x   5f         6x   6f
%
%  -3a -  2b -   c +   d = 1
%
%   9a +  4b +   c +   d = 0
%
% -27a -  8b -   c +   d = 0
%
%  81a + 16b +   c +   d = 0
%
%  Simultaneous solution for a, b, c and d followed by the solu-
%  tion of the preceding Taylor series, truncated after the u
%                                                            4x
%  terms, for u1  gives the following five-point approximation
%               x
%                                                                4
%  un-1  = (1/12dx)(-un-4 + 6un-3 - 18un-2 + 10un-1 + 3un) + O(dx )
%      x
%                                                              (4)
%
%    (5)  Right end, point i = n
%
%                                       2             3
%  a(un-4 = un + un (-4dx)  + un  (-4dx)  + un  (-4dx)
%                  x   1f       2x   2f       3x   3f
%
%                         4             5             6
%             + un  (-4dx)  + un  (-4dx)  + un  (-4dx)  + ...)
%                 4x   4f       5x   5f       6x   6f
%
%                                       2             3
%  b(un-3 = un + un (-3dx)  + un  (-3dx)  + un  (-3dx)
%                  x   1f       2x   2f       3x   3f
%
%                         4             5             6
%             + un  (-3dx)  + un  (-3dx)  + un  (-3dx)  + ...)
%                 4x   4f       5x   5f       6x   6f
%
%                                       2             3
%  c(un-2 = un + un (-2dx)  + un  (-2dx)  + un  (-2dx)
%                  x   1f       2x   2f       3x   3f
%
%                         4             5             6
%             + un  (-2dx)  + un  (-2dx)  + un  (-2dx)  + ...)
%                 4x   4f       5x   5f       6x   6f
%
%                                       2             3
%  d(un-1 = un + un ( -dx)  + un  ( -dx)  + un  ( -dx)
%                  x   1f       2x   2f       3x   3f
%
%                         4             5             6
%             + un  ( -dx)  + un  ( -dx)  + un  ( -dx)  + ...)
%                 4x   4f       5x   5f       6x   6f
%
%   -4a -  3b -  2c -   d = 1
%
%   16a +  9b +  4c +   d = 0
%
%  -64a - 27b -  8c -   d = 0
%
%  256a + 81b + 16c +   d = 0
%
%  Simultaneous solution for a, b, c and d followed by the solu-
%  tion of the preceding Taylor series, truncated after the u
%                                                            4x
%  terms, for u1  gives the following five-point approximation
%               x
%                                                                4
%  un  = (1/12dx)(3un-4 - 16un-3 + 36un-2 - 48un-1 + 25un) + O(dx )
%    x
%                                                              (5)
%
%  The weighting coefficients for equations (1) to (5) can be
%  summarized as
%
%             -25   48  -36   16   -3
%
%              -3  -10   18   -6    1
%
%       1/12    1   -8    0    8   -1
%
%              -1    6  -18   10    3
%
%               3  -16   36  -48   25
%
%  which are the coefficients reported by Bickley for n = 4, m =
%  1, p = 0, 1, 2, 3, 4 (Bickley, W. G., Formulae for Numerical
%  Differentiation, Math. Gaz., vol. 25, 1941.  Note - the Bickley
%  coefficients have been divided by a common factor of two).
%
%  Equations (1) to (5) can now be programmed to generate the
%  derivative u (x) of function u(x).
%              x
%
%  Compute the spatial increment
   dx=(xu-xl)/(n-1);
   r4fdx=1./(12.*dx)
   nm2=n-2;
%
%  Equation (1) (note - the rhs of equations (1), (2), (3), (4)
%  and (5) have been formatted so that the numerical weighting
%  coefficients can be more easily associated with the Bickley
%  matrix above)
   ux(  1)=r4fdx*...
     ( -25.*u(  1) +48.*u(  2) -36.*u(  3) +16.*u(  4)  -3.*u(  5));
%
%  Equation (2)
   ux(  2)=r4fdx*...
     (  -3.*u(  1) -10.*u(  2) +18.*u(  3)  -6.*u(  4)  +1.*u(  5));
%
%  Equation (3)
   for i=3:nm2
     ux(  i)=r4fdx*...
     (  +1.*u(i-2)  -8.*u(i-1)  +0.*u(  i)  +8.*u(i+1)  -1.*u(i+2));
   end
%
%  Equation (4)
   ux(n-1)=r4fdx*...
     (  -1.*u(n-4)  +6.*u(n-3) -18.*u(n-2) +10.*u(n-1)  +3.*u(  n));
%
%  Equation (5)
   ux(  n)=r4fdx*...
     (   3.*u(n-4) -16.*u(n-3) +36.*u(n-2) -48.*u(n-1) +25.*u(  n));
  end

Answer 1

[Torsten]

In the script part:

global nr nz dr dz drs dzs r z Dc Dt ca Tk cae Tke h k E R rk0 v rho Cp Tw dH ncall

In pde_2:

global nr nz dr dz drs dzs r r0 z zl Dc Dt ca Tk cae Tke h k E R rk0 v rho Cp Tw dH ncall

You can see that r0 is not set. This makes dx = empty in dss004.

[Shahid Hasnain]

@Torsten Many thanks