How can i solve a system coupled of PDE with an ODE using finite difference?

I tried several methods, but i couldn't find the solution. for instance, i used Crunk-Nicelson finite difference method like following script but i don't know how can i apply the secend eq.(ODE) inside the matrix. the big problem here is that each incerment is too small.

% The Crank Nicolson Method
clc, clear, close
eps=0.5; tor=3; lf=5*10^-6;%microM    dm=60*8.6*10^-6; %m2/s    kg= 1.87*60; % m/s   area=0.6 ;% m^-1
cin=10; %   ki=1*60; %mg/m3s  Ki=0.24; %m3/mg  miu=0.3*10^-6; %microM^-1  beta=0.5;   iif=16.5*10;%mw/cm2 
Kw=4.9*10^-4;% m^3/mg    cw=1000; %mg/m3  ux=0.011*60; %m/s
% syms x   % d=int((exp(-miu*lf*x)^beta)/lf,0,lf);
%%%%%%%%%%%%%%%%%%parameter
f=iif^beta; g=ki*Ki; h=Kw*cw; alfa=eps*dm/tor;
tf=180; nx=100;   dx=lf/nx; nt=2000;   dt=tf/nt;
% --- Constant Coefficients of the tridiagonal system
c = alfa/(2*dx^2)+ux/(4*dx);  % Subdiagonal: coefficients of u(i-1)
b = alfa/(2*dx^2)-ux/(4*dx);  % Super diagonal: coefficients of u(i+1)
a = 1/dt+b+c;  % Main Diagonal: coefficients of u(i)
% Boundary conditions and Initial Conditions
Uo(1)=20; Uo(2:nx)=0; 
Un(1)=20; Un(nx)=0;
% Store results for future use
UUU(1,:)=Uo;
% Loop over time
for k=2:nt
for ii=1:nx-2
if ii==1
d(ii)=c*Uo(ii)+(1/dt-b-c)*Uo(ii+1)+b*Uo(ii+2)+c*Un(1);%-f*g*Uo(ii+1)/(1+h+Ki*Uo(ii+1));
elseif ii==nx-2
d(ii)=c*Uo(ii)+(1/dt-b-c)*Uo(ii+1)+b*Uo(ii+2)+b*Un(nx);%-f*g*Uo(ii+1)/(1+h+Ki*Uo(ii+1));
else
d(ii)=c*Uo(ii)+(1/dt-b-c)*Uo(ii+1)+b*Uo(ii+2);%-f*g*Uo(ii+1)/(1+h+Ki*Uo(ii+1)); 
end
end  % note that d is row vector
%%%%%%%%%%%%%%%%%
% Transform a, b, c constants in column vectors:
bb=b*ones(nx-3,1);
cc=c*ones(nx-3,1);
aa=a*ones(nx-2,1);
% Use column vectors to construct tridiagonal matrices
AA=diag(aa)+ diag(-bb,1)+ diag(-cc,-1);
% Find the solution for interior nodes i=2,3,4,5
% UU=AA\d';
UU=inv(AA)*d';
% Build the whole solution as row vector
Un=[Un(1),UU',Un(nx)];
Un(nx)=Un(nx-1);
UUU(k,:)=Un;
Uo=Un; 
end
t=[0:dt:tf-dt];
uf=UUU(:,end);
plot(t,uf)

[Kuldeep Malik]

What if initial and boundary conditions are functions instead of constants

Answer 1

[Torsten]

Take a look at the answer provided here:

https://de.mathworks.com/matlabcentral/answers/371313-error-in-solving-system-of-two-reaction-diffusion-equations

The problem is very similar to yours.

Best wishes

Torsten.

[Moji]

i tried to apply that code for my probblem but i couldn't reach the answer. can you please find my mistake?

[Torsten]

You set

dC_iB/dx = 0

as boundary condition at x=0 instead of

C_iB = C_iB0

Best wishes

Torsten.

[Moji]

Dear Torsten I want to know how can i discretize CB and Cs in space that can be coupled. actualy, i couldn't undrestand that code you refrenced me. Thanks in advance

[Torsten]

Replace

u = zeros(n,1);
u(1,1)=20;

by

u = zeros(n,1);
u(1,1)=20;
u(n,1)=?

Here, 20 and ? should be the values for Ci_b at x=0 and x=L_f.

Then, in "fun", set

DuDt(1) = 0.0

before the line

DuDt(n)=0.

But setting values for Ci_b at inlet and outlet doesn't seem correct in my opinion. Further, I didn't check how you included the integral expression.

If the problem you posted is correct, the only difficulty is to approximate first and second order spatial derivatives in the inner grid points for Ci_b. This is done in the loop

for i=2:n-1
  DuDx(i) = ((x(i)-x(i-1))/(x(i+1)-x(i))*(u(i+1)-u(i))+(x(i+1)-x(i))/(x(i)-x(i-1))*(u(i)-u(i-1)))/(x(i+1)-x(i-1));
  D2uDx2(i) = (xhalf(i)*(u(i+1)-u(i))/(x(i+1)-x(i))-xhalf(i-1)*(u(i)-u(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));
end

With the help of these derivatives, it should be no problem to constuct the time derivatives DuDt and DvDt of Ci_B and Ci_S in the grid points.

Best wishes

Torsten.

[Moji]

B.C at the end point(nx+1) is unknown, i know u(nx+1)=u(nx). in addition for second dependent variable(Cs), i just have initial value. actually, i took this loop function for DuDx(i) and D2uDx2(i) from another code that i dont know how it is reached. my problem exactly is here how should i define DuDx(i) and D2uDx2(i)?? it am dizzy now.

[Torsten]

Sorry, I don't understand your problem.

The loop as it stands already calculates approximations to first and second order spatial derivatives for the function u (which is Cb_i in your case) in the inner grid points. No need to "define" DuDx(i) and D2uDx2(i).

I think it's time that you consult a textbook on the "method-of-lines" approach to solve partial differential equations.

Best wishes

Torsten.

[Moji]

Dear Mr. Torsten, i could run the problem without any error. But the answer is not correct (as i uploud it in attached file). i attached my final code here, could you please run it and tell me what is my mistake. the problem here is that solution gets the final answer very fast and time doesn't have any effect on answer.

Thanks a lot. Mojtaba

[Torsten]

The equations from above still hold ? Because you wrote that the boundary condition at the end point is unknown whereas in your .JPG file, the boundary values are known and fixed.

Best wishes

Torsten.

[Torsten]

Your new code is too much different from the one I linked to.

Sorry, but I don't have the time to dive in that deep.

Best wishes

Torsten.

[Torsten]

As you can see from the Username, it's code I wrote by myself.

Maybe if you ask more clearly what you don't understand, I will be able to explain.

Best wishes

Torsten.