ode45 concentration problem

Hello friends I wrote this code because,. But when I run the code, plot is wrong

Information about the system:

We have a lake which have zero concentration of organic pollutant C0=0mg/L, but there is an inflow to the lake which contain organic pollutant with a concentration Cin=31mg/L , when we measure the concentration of organic pollutant in the lake when t=15day=5475hr

C5475=12mg/L. I guess It should increase to the Cin=31mg/L

Cin= Concentration which is inflow of the system(lake), it is consumed with a rate k=5*(10^-6)

C=concentration of the lake, because of being complete mix system C=Cout.

V=volume of the lake

Q=inflow and ourflow of the lake

Min=total mass flux which also equal to Q*Cin

r=k.Cin*exp(-k*t)^2 which is second order decay reaction term

Actual complete mix system formula:

V*(dC/dt)=((Q*Cin)-(Q*C)-(r*V))

How to write a proper code for my purpose? Could you help me please, I need your help?

clear all;
clc;
close all;
%15day=365hr*15day=5475hr

C0=0;
C5475=12;
tspan = [0 5475]; 
[t,C]=ode45(@concentration, tspan, C0);

plot(t,C)
xlabel('time (hr)')
ylabel('Concentration (mg/L)')

function dCdt=concentration(t,C)
k=5*(10^-6);
A=100*10000;
h=2.5;
V=(A*h)/1000;
Q=500000;
Min=15.5*1000000;
Cin=Min/Q;
dCdt=((Q*Cin)/V-(Q*C)/V-(k*(Cin*exp(-k*t)^2)));
end

[Ethan Boyce]

In your function, I believe your concentration should be time dependent. I am not familiar with the ODE specifically but at first glance it seems like you're missing the t variable somewhere

[Carey n'eville]

I fixed it thank you

[Cris LaPierre]

I would think we'd need more information about the system. At the least some details about rates in and out, excretion/consumption, etc. Is the model based on a known system (e.g. three compartment model)?

[Carey n'eville]

I completed informations

[Cris LaPierre]

Units? Specifically

• V
• Q
• k

Answer 1

[Cris LaPierre]

First big error I see is that, assuming your A and h are in meters, you have not correctly converted your to L. . You have divided where you should have multiplied. Try this.

C0=0;
tspan = [0 5475]; 
[t,C]=ode45(@concentration, tspan, C0);

plot(t,C)
xlabel('time (hr)')
ylabel('Concentration (mg/L)')

function dCdt=concentration(t,C)
k=5E-6;
A=100*10000;
h=2.5;
V=(A*h)*1000; % Changed to multiply
Q=500000;
Cin=31; % mg/L
dCdt=((Q*Cin)/V-(Q*C)/V-(k*(Cin*exp(-k*t)^2)));
end

[Carey n'eville]

I appreciate to you your answer, it works thank you so much!