# MATLAB steady state (bvp4c) and transient (pdepe) plots do not match

I'm trying to use MATLAB solve a reaction-diffusion equation in both steady state (using bvp4c) and time dependency (using pdepe). When plotting both of the solutions, I do not get the same answer. For my boundary conditions I have: at x = 0, C = 0 and at x = L, dC/dx = 0.

For my steady state case (bvp4c):

  function steady_state_bvp4c
close all; clear all; clc;

%Diffusion coefficient
D_ij = 30; %3.0*10^-7 cm^2/s -> 30 um^2/s

%Initial concentration at x=0
L0 = 50;

%Total length
x_f = 20;

%Rate constants
k_1 = 0.00193;
k_2 = 0.00255;
k_3 = 4.09;
d_1 = 0.007;
d_2 = 3.95*10^-5;
d_3 = 2.26;

%Equilibrium constants
K_1 = 3.63;
K_2 = 0.0155;
K_3 = 0.553;
K_4 = 9.01;
K_i = 0.139;

%Total Receptors
N_T =1.7;

solinit = bvpinit(linspace(0,x_f,11),[0.5 0]);

sol = bvp4c(@(x,y)odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T),@twobc,solinit);

figure(1)
plot(sol.x,(sol.y(1,:)),'LineWidth',1)
xlabel('Distance from Device (\mum)')
ylabel('Concentration (nM)')
axis([0 x_f 0 L0])

function dy = odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T)
C_L    = y(1);
dC_Ldx = y(2);

N_r = -(K_3.*K_4.*K_i.*(K_2 + C_L - ((8.*K_2.*C_L.^2.*N_T + 8.*K_2.*K_3.*C_L.*N_T + K_2.^2.*K_3.*K_4.*K_i + K_3.*K_4.*K_i.*C_L.^2 + 8.*K_2.^2.*K_3.*K_i.*N_T + 2.*K_2.*K_3.*K_4.*K_i.*C_L + 8.*K_2.*K_3.*K_i.*C_L.*N_T)./(K_3.*K_4.*K_i)).^(1./2)))./(4.*C_L.^2 + 4.*K_3.*C_L + 4.*K_2.*K_3.*K_i + 4.*K_3.*K_i.*C_L);
N_R = (C_L./K_1).*N_r;
N_r2 = (1./K_4).*(N_r).^2;
N_rR = (C_L./(2.*K_2.*K_4)).*(N_r).^2;
N_pR = (C_L./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
N_R2 = ((C_L.^2)./(K_2.*K_3.*K_4.*K_i)).*(N_r).^2;

R_L = ((2.*d_1.*(N_T-N_r-N_r2-N_rR-N_pR-N_R2))-(2.*k_1.*C_L.*N_r))+...
((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*C_L.*(N_T-N_r-N_R-N_rR-N_pR-N_R2)))+...
(d_3.*(N_T-N_r-N_R-N_r2-N_rR-N_pR))-(2.*k_3.*C_L.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));

dy = zeros(2,1);
dy(1) = dC_Ldx;
dy(2) = (-R_L/D_ij);
end

function res = twobc(ya,yb)
res = [ ya(1)-(L0); yb(2) ];
end
end


For my time dependent case:

  function time_dependent_pdepe
clear all; close all; clc;

D_ij = 30; %Diffusion coefficient D (3.0*10^-7 cm^2/s -> 30 um^2/s)
L0 = 50; %c0 [nM]
x_f =20; %Length of domain [um]
maxt = 1; %Max simulation time [s]

m = 0; %Parameter corresponding to the symmetry of the problem
x = linspace(0,x_f,100); %xmesh
t = linspace(0,maxt,100); %tspan

sol = pdepe(m,@DiffusionPDEfun,@DiffusionICfun,@DiffusionBCfun,x,t,[]);
u = sol;

% Plotting
hold all
for n = linspace(1,length(t),10)
plot(x,sol(n,:),'LineWidth',2)

end
title('Time Dependent')
xlabel('Distance from Device (\mum)')
ylabel('Concentration (nM)')
axis([0 x_f 0 L0])

function [c,f,s] = DiffusionPDEfun(x,t,u,dudx)

D = D_ij;

%Rate constants
k_1 = 0.00193;
k_2 = 0.00255;
k_3 = 4.09;
d_1 = 0.007;
d_2 = 3.95;
d_3 = 2.26;

%Equilibrium constants
K_1 = 3.63;
K_2 = 0.0155;
K_3 = 0.553;
K_4 = 9.01;
K_i = 0.139;

%Total Recptors
N_T = 1.7;

N_r = -(K_3.*K_4.*K_i.*(K_2 + u - ((8.*K_2.*u.^2.*N_T + 8.*K_2.*K_3.*u.*N_T + K_2.^2.*K_3.*K_4.*K_i + K_3.*K_4.*K_i.*u.^2 + 8.*K_2.^2.*K_3.*K_i.*N_T + 2.*K_2.*K_3.*K_4.*K_i.*u + 8.*K_2.*K_3.*K_i.*u.*N_T)./(K_3.*K_4.*K_i)).^(1./2)))./(4.*u.^2 + 4.*K_3.*u + 4.*K_2.*K_3.*K_i + 4.*K_3.*K_i.*u);
N_R = (u./K_1).*N_r;
N_r2 = (1./K_4).*(N_r).^2;
N_rR = (u./(2.*K_2.*K_4)).*(N_r).^2;
N_pR = (u./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
N_R2 = ((u.^2)./(K_2.*K_3.*K_4.*K_i)).*(N_r).^2;

R_L = ((2.*d_1.*(N_T-N_r-N_r2-N_rR-N_pR-N_R2))-(2.*k_1.*u.*N_r))+...
((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*u.*(N_T-N_r-N_R-N_rR-N_pR-N_R2)))+...
(d_3.*(N_T-N_r-N_R-N_r2-N_rR-N_pR))-(2.*k_3.*u.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));

c = 1;
f = D_ij.*dudx;
s = R_L;
end

function u0 = DiffusionICfun(x)
u0 = 0;
end

function [pl,ql,pr,qr] = DiffusionBCfun(xl,ul,xr,ur,t)
% Boundary conditions for x = 0 and x = L;

c0 = L0;

% BCs: No flux boundary at the right boundary and constant concentration on
% the left boundary

% At x = 0, c0 = L0
pl = ul-c0;
ql = 0;

% At x = L, dc/dx = 0
pr = 0;
qr = 1;

end

end


time dependent plot

Is there a difference in the way I am setting up both situations? I would greatly appreciate any help! Thank you!

• There is way to much code for someone to go through, so you are unlikely to get an answer. I would suggest finding the absolute simplest setting where you get different answers for time-dependent and stationary code, and post that. Doing this type of exercise will probably illuminate the issue for you anyway. – Jeff Jul 17 '18 at 17:43