function time_dependent_pdepe_NEW
clear all; close all; clc;
D_ij = 30;
L0 = 100;
x_f =1000;
maxt = 1000;
m = 0;
x = linspace(0,x_f,100);
t = linspace(0,maxt,10);
sol = pdepe(m,@DiffusionPDEfun,@DiffusionICfun,@DiffusionBCfun,x,t,[]);
u = sol;
hold all
for n = linspace(1,length(t),10)
plot(x,sol(n,:),'LineWidth',2)
end
function [c,f,s] = DiffusionPDEfun(x,t,u,dudx)
D = D_ij;
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;
K_1 = 3.63;
K_2 = 0.0155;
K_3 = 0.553;
K_4 = 9.01;
K_i = 0.139;
N_T = 1.7;
N_r = (-1 - (u./K_1) + sqrt((-1 - (u./K_1)).^2 -4.*(-(2./K_4) - ((2.*u)./(K_2.*K_4)) - ((2.*u)./(K_2.*K_4.*K_i)) - ((2.*u.^2)/(K_2.*K_3.*K_4.*K_i))).*N_T))/(4.*((1./K_4) + (u./(K_2.*K_4)) + (u./(K_2.*K_4.*K_i)) + (u.^2./(K_2.*K_3.*K_4.*K_i))));
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 = N_T-(2.*N_R2)-(4.*N_pR)-N_r-(4.*N_rR)-(2.*N_r2);
r_2 = (N_T./2)-(N_R./2)-N_R2-(2.*N_pR)-(N_r./2)-(2.*N_rR);
rR = (N_T./4)-(N_R./4)-(N_R2./2)-N_pR-(N_r./4)-(N_r2./2);
pR = (N_T./4)-(N_R./4)-(N_R2./2)-(N_r./4)-N_rR-(N_r2./2);
R_2 = (N_T./2)-(N_R./2)-(2.*N_pR)-(N_r./2)-(2.*N_rR)-N_r2;
R_L = ((2.*d_1.*(R))-(2.*k_1.*u.*N_r))+...
((2.*d_2.*(rR)))-((k_2.*u.*(r_2)))+...
(d_3.*(R_2))-(2.*k_3.*u.*(pR));
c = 1;
f = D_ij.*dudx;
s = R_L;
end
function u0 = DiffusionICfun(x)
c0 = L0;
if(x==0)
u0 = c0;
else
u0 = 0;
end
end
function [pl,ql,pr,qr] = DiffusionBCfun(xl,ul,xr,ur,t)
c0 = L0;
pl = ul-c0;
ql = 0;
pr = 0;
qr = 1;
end
end
