Comparing two row vectors to find constant slope for steady state condition

I am comparing two rows to find if at any point I am getting a uniform slope (a condition for steady state let's say in a diffusion problem)

c2(N+1,N+1)=9.90083146790674e-07;    % Doing this to add an extra index so I can loop through this 100*100 matrix 

for i = 1:N 
    p(:,i) = c2(:,i)-c2(:,i+1); 
    if p(:,i) < 1e-10         % I wany to see if the slope is uniform/0 so I can ascertain if my solution has reached steady state (for a diffusion problem) 
        disp(i)
    end
end

Thanks!

[Mathieu NOE]

hello

so basically you are doing a difference (like diff) and compare that to a threshold.

That's a valid approach to find a constant slope condition as you say (steady state)

so what is the issue ? what are you willing to do ?

[Hashim]

Hi Sorry I got the wrong variable here... it's actually c2 (see attached file) that I am concerned with. if you run the code you will find it output a graph (figure 2) with two y axis. I wish to know at what time the slope on the right curve is removed or lessened signifncatly so I can know that my system is at a steady state.

Answer 1

[Mathieu NOE]

Hi again

so I introduced the first and second derivative of c2 and plotted in figure 3

Attached is the function to do a finite difference derivation better than using diff (or alike)

as the second derivative goes to zero in the second half of the x data, you can tell the system has reached a steady state

main code (modified) below :

% function pdepe_philip_v7a

% 5- We now add a second order reaction to the mix and try to study the
% effects of this reaction on the concentration output of the system 
% Also, to understand the concept of reaction layers and how they can explain 
% behavior or our system 

% Solve for Substrate concentration as well 

global D_M D_S M_bulk n F m Area R S_bulk k1

n          = 1;                                  % No. of Electrons Involved
F          = 96485.3329;                         % sA/mol
R          = 8.3145;                             % kgcm^2/s^2.mol.K
D_M        = 5e-06;                              % cm^2/s 
D_S        = 5e-06;                              % cm^2/s (???)
l          = 0.01;                               % cm 
Area       = 1;                                  % cm^2 
M_bulk     = 1e-06;                              % mol/cm^3
S_bulk     = 1e-04;                              % mol/cm^3 

m          = 0;                                  % Cartesian Co-ordinates

k1         = 1e+05; 

N          = 100; 
 
% computational cost of the solution depends weakly on the length of tspan
t          = linspace(0, 10, N);                % s       

% xmesh (Determines the computational cost)
x          = linspace(0, l, N);

sol        = pdepe(m, @pdev7, @pdev7ic, @pdev7bc, x, t);
c1         = sol(:, :, 1);                        % Mox Conc.  
c2         = sol(:, :, 2);                        % Substrate Conc.
c3         = sol(:, :, 3);                        % Mred Conc.  

% % Conc. Profiles (3D)
% figure(1)
% surf(x,t,c2);       
% view(0, 0);
% xlabel('x');
% ylabel('t(sec)');
% zlabel('Concentration');

% % dc/dt [Glucose] 
% figure(3); 
% plot(x, c1( 2:N, :));  




SS         = 94;

% dc/dt [Os] + [Glucose]
figure(2); 
plotyy(x, c1( SS, :), x, c2( SS,: )); 
xlabel('x (cm)'); ylabel('[Os^3] (mol/cm^3)'); 
hold on;  

% draw a line for sigma_kinetic/reaction layer [cm]
x_k        = sqrt(((D_S*M_bulk)+(D_M*S_bulk))/...
             (k1*(M_bulk+S_bulk).^2));     
         
if S_bulk > M_bulk
    xl          = [x_k, x_k];
    yl          = [-M_bulk*2, M_bulk*2];
    line(xl, yl, 'Color','green','LineStyle','--','LineWidth',2);
else 
    xl          = [l-x_k, l-x_k];
    yl          = [-M_bulk*2, M_bulk*2];
    line(xl, yl, 'Color','g','LineStyle','--','LineWidth',2);
end
hold off;  

%%  plot of first and second derivative of C2
  [dc2, ddc2] = firstsecondderivatives(x,c2( SS,: ));
figure(3); 
subplot(3,1,1),plot(x,c2( SS,: ))
xlabel('x (cm)'); ylabel('C2'); 
subplot(3,1,2),plot(x,dc2)
xlabel('x (cm)'); ylabel('dC2/dx'); 
subplot(3,1,3),plot(x,ddc2)
xlabel('x (cm)'); ylabel('d²C2/dx²'); 

% % Mred Flux at x=0
% figure(4);
% subplot(4,1,1);
% plot(t, c1( :,1)); 
% xlabel('time (s)'); ylabel('dM/dx @ x=0)');  
% 
% % Mred Flux at x=d
% subplot(4,1,2);
% plot(t, c1(:,N));
% xlabel('time (s)'); ylabel('dM/dx @ x=d)');  
% 
% % Substrate Flux at x=d
% subplot(4,1,3);
% plot(t, c2(:,N));
% xlabel('time (s)'); ylabel('dS/dx @ x=d)');  
% 
% % Mred Flux at x=0 - Mred Flux at x=0
% subplot(4,1,4);
% plot(t, c1(:,N)-c1(:,1));
% xlabel('time (s)'); ylabel('dM/dx@x=d-dM/dx@x=0');  

asd  = c1(:,N-1) - c1(:,N); 

c1(N+1,N+1)=9.90083146790674e-07;

for i = 1:N 
    p(:,i) = c1(:,i)-c1(:,i+1); 
    if p(:,i) < 1e-10
        disp(i)
    end
end

as = 1

% end

%% pdepe Function That Calls Children
%%

function [a, f, s] = pdev7(x, t, c, DuDx)

global D_M D_S k1  
a          = [1; 1; 1]; 
f          = [D_M; D_S; D_M].*DuDx;    

% c(1)--> Mox(x,t) || c(2)--> S(x,t) || c(3)--> Mred(x,t)
s          = [-(k1*c(1)*c(2)); -k1*c(1)*c(2); k1*c(1)*c(2)] ;

end

%% Initial Condition
%%

function c0 = pdev7ic(x)

global S_bulk M_bulk 
c0          = [M_bulk; S_bulk; M_bulk];
end

%% Boundry Conditions 
%%

function [pl, ql, pr, qr] = pdev7bc(xl, cl, xr, cr, t)

global M_bulk S_bulk n F R
E0   =0.2; 
E    =1.2; 
alpha=exp((E-E0)*((n*F)/(R*298.15)));   
pl   =[cl(1)-((M_bulk*alpha)./(1+alpha)); 0; cl(3)-(M_bulk*(1/(1+alpha)))];
ql   =[0;                                    1;                            0];
pr   =[0;                                    cr(2)-S_bulk;                 0];
qr   =[1;                                    0;                            1];   

end



%% Analytical Portion Solution 
%%

[Hashim]

Hi Mathieu,

I am really really grateful that you took the time to go through my question.

So, as per my understanding what you've done is you're taking a derivate of c2 which, and please correct me if I am wrong, is the ( 𝜕𝑢)/𝜕𝑥 from the standard form 𝑔(𝑥, 𝑡, 𝑢 ( 𝜕𝑢)/𝜕𝑥) 𝜕𝑢/𝜕𝑡=𝑥^(−𝑚) 𝜕/𝜕𝑥 (𝑥^𝑚 𝑓(𝑥, 𝑡, 𝑢 ( 𝜕𝑢)/𝜕𝑥)) as understood by pdepe. And, if it is then isn't your ( 𝜕2𝑢)/𝜕𝑥2 really ( 𝜕3𝑢)/𝜕𝑥3 since you've already taken derivate two times before.

[Mathieu NOE]

hello

I simply here considered c2 as function of x sollely and derived c2 once and twice vs x

as the second derivative becomes small => steady state is reached

hope it helped

all the best

[Hashim]

Would it be alright if you allow me a few days so I can consult my advisor on this before I go on to accept your answer. I have a few presentations in a week that's why.

[Mathieu NOE]

Ok no problem !

[Hashim]

Hi @Mathieu NOEMathew I've got the code working but I need to see the variable against time and not x axis. I've done that in the code below. Now I would like to check if it works by taking a snapshot of system when it is at an earlier time. For this I need to plot t(some earlier value of time) against c2(some earlier value of time). Now we know that c2 is a square matrix and time is a vector. How do I code it?

I am talking about these line BTW in the code below.

% plot of first and second derivative of C2

[dy, ddy] = firstsecondderivatives(t, c2(SS,:));

figure(3);

subplot(3,1,1),plot(t,c2(SS,: ));

xlabel('t (s)'); ylabel('C2');

subplot(3,1,2),plot(t,dy)

xlabel('t (s)'); ylabel('dC2/dx');

subplot(3,1,3),plot(t,ddy)

xlabel('t (s)'); ylabel('d²C2/dx²');

function pdepe_philip_v7a

% 5- We now add a second order reaction to the mix and try to study the
% effects of this reaction on the concentration output of the system 
% Also, to understand the concept of reaction layers and how they can explain 
% behavior or our system 

% Solve for Substrate concentration as well 

global D_M D_S M_bulk n F m Area R S_bulk k1

n          = 1;                                  % No. of Electrons Involved
F          = 96485.3329;                         % sA/mol
R          = 8.3145;                             % kgcm^2/s^2.mol.K
D_M        = 5e-06;                              % cm^2/s 
D_S        = 5e-06;                              % cm^2/s (???)
l          = 0.01;                               % cm 
Area       = 1;                                  % cm^2 
M_bulk     = 1e-06;                              % mol/cm^3
S_bulk     = 1e-04;                              % mol/cm^3 

m          = 0;                                  % Cartesian Co-ordinates

k1         = 1e+05; 

N          = 100; 
 
% computational cost of the solution depends weakly on the length of tspan
t          = linspace(0, 100, N);                % s       

% xmesh (Determines the computational cost)
x          = linspace(0, l, N);

sol        = pdepe(m, @pdev7, @pdev7ic, @pdev7bc, x, t);
c1         = sol(:, :, 1);                       % Mox Conc.  
c2         = sol(:, :, 2);                       % Substrate Conc.
c3         = sol(:, :, 3);                       % Mred Conc.  

% % Conc. Profiles (3D)
% figure(1)
% surf(x,t,c2);       
% view(0, 0);
% xlabel('x');
% ylabel('t(sec)');
% zlabel('Concentration');

SS         = 10;

% dc/dt [Os] + [Glucose]
figure(2); 
[hAx,hLine1,hLine2]= plotyy(x, c1( SS, :), x, c2( SS,: )); 
xlabel('x / cm'); 
ylabel(hAx(1),'Os^3^+ / mol cm^-^3')    % left y-axis 
ylabel(hAx(2),'Substrate / mol cm^-^3') % right y-axis
% legend('k=1e5 cm^3 mol s^-^1','D=5e-6 cm^2 s^-^1','t=100 s');
hold on;  

% draw a line for sigma_kinetic/reaction layer [cm]
x_k            = sqrt(((D_S*M_bulk)+(D_M*S_bulk))/...
                 (k1*(M_bulk+S_bulk).^2));     
         
if S_bulk > M_bulk
    x_new      = [x_k, x_k];
    y          = [-M_bulk*3, M_bulk*2];
    line(x_new, y, 'Color','green','LineStyle','--','LineWidth',2);
else 
    x_newer    = [l-x_k, l-x_k];
    y          = [-M_bulk*3, M_bulk*2];
    line(x_newer, y, 'Color','g','LineStyle','--','LineWidth',2);
end

size(c2(SS,1)) 
size(t)

% plot of first and second derivative of C2
[dy, ddy] = firstsecondderivatives(t, c2(SS,:));
figure(3); 
subplot(3,1,1),plot(t,c2(SS,: ));
xlabel('t (s)'); ylabel('C2'); 
subplot(3,1,2),plot(t,dy)
xlabel('t (s)'); ylabel('dC2/dx'); 
subplot(3,1,3),plot(t,ddy)
xlabel('t (s)'); ylabel('d²C2/dx²'); 

% % dMdx at x=0
% figure(4);
% subplot(4,1,1);
% plot(t, c1( :,1)); 
% xlabel('time (s)'); ylabel('dM/dx @ x=0)');  
% 
% % dMdx at x=d
% subplot(4,1,2);
% plot(t, c1(:,N));
% xlabel('time (s)'); ylabel('dM/dx @ x=d)');  
% 
% % dSdx at x=d
% subplot(4,1,3);
% plot(t, c2(:,N));
% xlabel('time (s)'); ylabel('dS/dx @ x=d)');  
% 
% % dMdx at x=0 - dMdx at x=0
% subplot(4,1,4);
% plot(t, c2(N,:)-c1(2,:));
% xlabel('time (s)'); ylabel('dM/dx@x=d-dM/dx@x=0');  

% asd  = c2(N-1,:) - c2(N,:); 
% 
% c2(N+1,N+1)=9.90083146790674e-07;
% 
% for i = 1:N 
%     p(i,:) = c2(i,:)-c2(i+1,:); 
%     if p(i,:) < 1e-8
%         disp(i)
%     end
% end

% slope   = gradient(c2(SS,:), linspace(0, l, N)); 


end

%% pdepe Function
%%

function [a, f, s] = pdev7(x, t, c, DuDx)

global D_M D_S k1  
a          = [1; 1; 1]; 
f          = [D_M; D_S; D_M].*DuDx;    

% c(1)--> Mox(x,t) || c(2)--> S(x,t) || c(3)--> Mred(x,t)
s          = [-(k1*c(1)*c(2)); -k1*c(1)*c(2); k1*c(1)*c(2)] ;

end

%% Initial Condition
%%

function c0 = pdev7ic(x)

global S_bulk M_bulk 
c0          = [M_bulk; S_bulk; M_bulk];
end

%% Boundry Conditions 
%%

function [pl, ql, pr, qr] = pdev7bc(xl, cl, xr, cr, t)

global M_bulk S_bulk n F R
E0   =0.2; 
E    =1.2; 
alpha=exp((E-E0)*((n*F)/(R*298.15)));   
pl   =[cl(1)-((M_bulk*alpha)./(1+alpha)); 0; cl(3)-(M_bulk*(1/(1+alpha)))];
ql   =[0;                                    1;                            0];
pr   =[0;                                    cr(2)-S_bulk;                 0];
qr   =[1;                                    0;                            1];   

end

%% Courtesy of Mathieu From MATLAB Website 
%%

function [dy, ddy] = firstsecondderivatives(x,y)

% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.
% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.

n      = length(x);
dy     = zeros;
ddy    = zeros;

% Input variables:
% x: vector with the x the data points.
% y: vector with the f(x) data points.
% Output variable:
% dy: Vector with first derivative at each point.
% ddy: Vector with second derivative at each point.
dy(1)  = (-3*y(1) + 4*y(2) - y(3)) / (2*(x(2) - x(1)));       % First derivative
ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative

for i  = 2:n-1   
	dy(i)  = (y(i+1) - y(i-1)) / (x(i+1) - x(i-1));
	ddy(i) = (y(i-1) - 2*y(i) + y(i+1)) / (x(i-1) - x(i))^2;
end
	dy(n)  = (y(n-2) - 4*y(n-1) + 3*y(n)) / (2*(x(n) - x(n-1)));
	ddy(n) = (-y(n-3) + 4*y(n-2) - 5*y(n-1) + 2*y(n)) / (x(n) - x(n-1))^2;
end

%%
%%

[Hashim]

Hi @Mathieu NOE I've been playing around with your function and pdeval which also provided dc/dx and found that even though the gradient is the same the values are off by a magnitude of 10s. Could you shine some light on this?

I've added another subplot where you can see this...

function pdepe_philip_v7a

% 5- We now add a second order reaction to the mix and try to study the
% effects of this reaction on the concentration output of the system 
% Also, to understand the concept of reaction layers and how they can explain 
% behavior or our system 

% Solve for Substrate concentration as well 

global D_M D_S M_bulk n F m Area R S_bulk k1

n          = 1;                                  % No. of Electrons Involved
F          = 96485.3329;                         % sA/mol
R          = 8.3145;                             % kgcm^2/s^2.mol.K
D_M        = 5e-06;                              % cm^2/s 
D_S        = 5e-06;                              % cm^2/s (???)
l          = 0.01;                               % cm 
Area       = 1;                                  % cm^2 
M_bulk     = 1e-06;                              % mol/cm^3
S_bulk     = 1e-04;                              % mol/cm^3 

m          = 0;                                  % Cartesian Co-ordinates

k1         = 1e+05; 

N          = 100; 
 
% computational cost of the solution depends weakly on the length of tspan
t          = linspace(0, 100, N);                % s       

% xmesh (Determines the computational cost)
x          = linspace(0, l, N);

sol        = pdepe(m, @pdev7, @pdev7ic, @pdev7bc, x, t);
c1         = sol(:, :, 1);                       % Mox Conc.  
c2         = sol(:, :, 2);                       % Substrate Conc.
c3         = sol(:, :, 3);                       % Mred Conc.  

SS         = 5; 

% dc/dt [Os] + [Glucose]
figure(2); 
[hAx,hLine1,hLine2]= plotyy(x, c1( SS, :), x, c2( SS,: )); 
xlabel('x / cm'); 
ylabel(hAx(1),'Os^3^+ / mol cm^-^3')    % left y-axis 
ylabel(hAx(2),'Substrate / mol cm^-^3') % right y-axis
% legend('k=1e5 cm^3 mol s^-^1','D=5e-6 cm^2 s^-^1','t=100 s');
 
% 
% % draw a line for sigma_kinetic/reaction layer [cm]
% x_k            = sqrt(((D_S*M_bulk)+(D_M*S_bulk))/...
%                  (k1*(M_bulk+S_bulk).^2));     
%          
% if S_bulk > M_bulk
%     x_new      = [x_k, x_k];
%     y          = [-M_bulk*3, M_bulk*2];
%     line(x_new, y, 'Color','green','LineStyle','--','LineWidth',2);
% else 
%     x_newer    = [l-x_k, l-x_k];
%     y          = [-M_bulk*3, M_bulk*2];
%     line(x_newer, y, 'Color','g','LineStyle','--','LineWidth',2);
% end

size(c2(1,1:SS)) 
size(t(1:SS))

% plot of first and second derivative of C2
% t1 = t(1:SS); 
c4 = c2(SS,:);
[dy, ddy] = firstsecondderivatives(t, c4);
figure(3); 
subplot(4,1,1),plot(t,c4);
xlabel('t (s)'); ylabel('C2'); 
subplot(4,1,2),plot(t,dy)
xlabel('t (s)'); ylabel('dC2/dx'); 
subplot(4,1,3),plot(t,ddy)
xlabel('t (s)'); ylabel('d²C2/dx²'); 
[cout,dcdx]= pdeval(m,x,c2(SS,:),x);
subplot(4,1,4),plot(t,dcdx)
xlabel('t (s)'); ylabel('dC2/dx');





% % dMdx at x=0
% figure(4);
% subplot(4,1,1);
% plot(t, c1( :,1)); 
% xlabel('time (s)'); ylabel('dM/dx @ x=0)');  
% 
% % dMdx at x=d
% subplot(4,1,2);
% plot(t, c1(:,N));
% xlabel('time (s)'); ylabel('dM/dx @ x=d)');  
% 
% % dSdx at x=d
% subplot(4,1,3);
% plot(t, c2(:,N));
% xlabel('time (s)'); ylabel('dS/dx @ x=d)');  
% 
% % dMdx at x=0 - dMdx at x=0
% subplot(4,1,4);
% plot(t, c2(N,:)-c1(2,:));
% xlabel('time (s)'); ylabel('dM/dx@x=d-dM/dx@x=0');  

% asd  = c2(N-1,:) - c2(N,:); 
% 
% c2(N+1,N+1)=9.90083146790674e-07;
% 
% for i = 1:N 
%     p(i,:) = c2(i,:)-c2(i+1,:); 
%     if p(i,:) < 1e-8
%         disp(i)
%     end
% end

% slope   = gradient(c2(SS,:), linspace(0, l, N)); 


end

%% pdepe Function
%%

function [a, f, s] = pdev7(x, t, c, DuDx)

global D_M D_S k1  
a          = [1; 1; 1]; 
f          = [D_M; D_S; D_M].*DuDx;    

% c(1)--> Mox(x,t) || c(2)--> S(x,t) || c(3)--> Mred(x,t)
s          = [-(k1*c(1)*c(2)); -k1*c(1)*c(2); k1*c(1)*c(2)] ;

end

%% Initial Condition
%%

function c0 = pdev7ic(x)

global S_bulk M_bulk 
c0          = [M_bulk; S_bulk; M_bulk];
end

%% Boundry Conditions 
%%

function [pl, ql, pr, qr] = pdev7bc(xl, cl, xr, cr, t)

global M_bulk S_bulk n F R
E0   =0.2; 
E    =1.2; 
alpha=exp((E-E0)*((n*F)/(R*298.15)));   
pl   =[cl(1)-((M_bulk*alpha)./(1+alpha)); 0; cl(3)-(M_bulk*(1/(1+alpha)))];
ql   =[0;                                    1;                            0];
pr   =[0;                                    cr(2)-S_bulk;                 0];
qr   =[1;                                    0;                            1];   

end

%% Courtesy of Mathieu From MATLAB Website 
%%

function [dy, ddy] = firstsecondderivatives(x,y)

% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.
% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.

n      = length(x);
dy     = zeros;
ddy    = zeros;

% Input variables:
% x: vector with the x the data points.
% y: vector with the f(x) data points.
% Output variable:
% dy: Vector with first derivative at each point.
% ddy: Vector with second derivative at each point.
dy(1)  = (-3*y(1) + 4*y(2) - y(3)) / (2*(x(2) - x(1)));       % First derivative
ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative

for i  = 2:n-1   
	dy(i)  = (y(i+1) - y(i-1)) / (x(i+1) - x(i-1));
	ddy(i) = (y(i-1) - 2*y(i) + y(i+1)) / (x(i-1) - x(i))^2;
end
	dy(n)  = (y(n-2) - 4*y(n-1) + 3*y(n)) / (2*(x(n) - x(n-1)));
	ddy(n) = (-y(n-3) + 4*y(n-2) - 5*y(n-1) + 2*y(n)) / (x(n) - x(n-1))^2;
end

%%
%%

[Mathieu NOE]

hello

I see that you do now derivative against t and not anymore against x . So I changed also the legend accordingly

the first and second derivatives against t are OK - see the rise of C , about 0.6 e- 5 over 100 seconds , so average slope is ~ 0.6 e-7.

if you compute the mean value of dy : mean(dy) = 6.4915e-08 , which is quite the same

now the 4th subplot is derivative against x , and delta x is not same as delta t (factor 10,000)

last but not least , my first derivative and gradient do give the same results : see simple example below :

x = 0:0.1:1;
y = 0:0.1:1;

gradient(x,y)  =      1     1     1     1     1     1     1     1     1     1     1

[dy, ddy] = firstsecondderivatives(x, y)

dy =

    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000


ddy =

   1.0e-13 *

         0         0         0   -0.0555         0         0         0    0.1110   -0.1110         0         0

[Mathieu NOE]

your figure (3) code updated :

figure(3); 
subplot(4,1,1),plot(t,c4);
xlabel('t (s)'); ylabel('C2'); 
subplot(4,1,2),plot(t,dy)
xlabel('t (s)'); ylabel('dC2/dt'); 
subplot(4,1,3),plot(t,ddy)
xlabel('t (s)'); ylabel('d²C2/dt²'); 
[cout,dcdx]= pdeval(m,x,c2(SS,:),x);
subplot(4,1,4),plot(t,dcdx)
xlabel('t (s)'); ylabel('dC2/dx');

[Hashim]

Ah yes sorry my bad that was quite obvious... I just chnged the t with x and the profiles match exactly. Would be even handy if pdeval gave the option to calcuate higher derivatives too. That being said what I've come to understand is that the pdepe function is giving me the conc. (c) values acorss the domain and not dc/dx as default. Also, the dc/dx tell me about the change in the conc. over the spatial domain... while on the other hand the second derivative d2c/dx2 tell me about the change in dc/dx.

Now in oder to ascertain when the system reaches the steady state I can use both dc/dx and d2c/dx2 but in different ways i.e.

dc/dx = constant slope means steady state

d2c/dx2 = 0 or something very small means steady state

Is that correct? Also how should one go about implementing it?

Grateful for your time.

[Mathieu NOE]

hello again

yes , basically "dc/dx = constant slope means steady state" is equivalent to " d2c/dx2 = 0 or something very small means steady state "

about implementing it : why not simply decide that if abs(d2c/dx2) < threshold (to be defined) over a long enough period (> 1s ) then steady state reached .

[Hashim]

How about this, anyways I've accepeted the answer, thank you extremely much for your help so far.

for i = 1:length(t) 
    if abs(ddy(i)) < 1e-9
        disp(i)
    end
end

[Mathieu NOE]

you can do the same without a for loop

ind = find(abs(ddy(i)) < 1e-9);
disp(t(ind))