Algae growing. Concentration curve problem

Hi! i have made the following code:

function w = Bioreactor

clc
clear all

% Forward Euler method is used for solving the following boundary value problem:

% Problem parameters
alpha1 =  5*10^(-4);                                       % Diffusion coefficient [m^2/s]
mu0 = 5; 
mue = 3.5; 
Hl = 5;
T = 5000;                                                  % [s]
nz = 100; nt = 1000;
KP = 2; 
HP= 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 

lambda = 2 ; 
rs= 10;

P(1)=3;                                                     % Initial chosen concentration of P [g/L]
N(1)=9;                                                     % Initial chosen concentration of N [g/L]
C(1)=30;                                                    % Initial chosen concentration of C [g/L]

Z = 10;                                                     % Lenght of the bioreactor [m]

%Check of the problem's parameters
if  any([P(1)-PC  N(1)-NC] <=0 ) 
    error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T  nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end

% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;

% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt);                                             % Round nt to an integer
nz = round(nz);                                             % Round nz to an integer
N = zeros(1, nt+1);                                         % Initialize N to zero
P = zeros(1, nt+1);                                         % Initialize P to zero
C = zeros(1, nt+1);                                         % Initialize C to zero

I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w = zeros(nz+1,nt+1);


% Finite-difference method
for j = 1:nt

    f1= (KP*(P(j)-PC))/(HP+(P(j)-PC));
    f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC));
    pH = (6.35 - log10(C(j))) /2 ;
    f3 = 1/(1+exp(lambda*(pH-PHs3)));
    f4 = 1/(1+exp(lambda*(pH-PHs4))); 
    
    I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
    g = mu0*I_in./(Hl+I_in);
    
    integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j));    % Data organization
    N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j));           % Function of N over the time and space
    P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j));           % Function of P over the time and space
    C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j));           % Function of C over the time and space
    
    w(nz+1,j+1) = 9;                                        % Algae concentration at z=0  [Dirichlet]
    w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + alpha1 * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2);
    w(1,j+1) = w(2,j+1);                                    % Algae concentration at z=10 [Neumann]
    
    
    plot(z,flipud(w(:,j+1)),'b');
    %set(gca,'XDir','reverse');
    xlabel('z'); ylabel('w');
    title(['t=',num2str(t(j)),' s']); 
    axis([0 10 0 10]);
    pause(0.0001);
end
end

And the code is working perfectly. But the W which is the algae concentration over space and time it has been fixed a value of 9 mg/L at the beginning of this bioreactor, so at z=0, im talking about the following row:

w(nz+1,j+1) = 9;                                        % Algae concentration at z=0  [Dirichlet]

But now i would like to fix this concentration at a certain z which has to be equal to z=Z/2 so at the half of the bioreactor. I know that in order to make the code working for sure i have to modify this row:

    w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + alpha1 * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2);

Can you help me to do that? I've tried something but anything really worked.

Thanks, attached to this post there is the file from which i took the equations.

I should have a graph like this:

[Alfonso]

@Torsten

Answer 1

[Torsten]

Here is a code for a simple heat-conduction equation

dT/dt = d^2T/dz^2
T(z,0) = 0
T(0,t) = 200
T(1,t) = 340
T(0.5,t) = 90

You should be able to adapt it to your application.

L = 1;
zcut = L/2;
n1 = 50;
n2 = 50;
z1 = linspace(0,zcut,n1);
dz1 = z1(2)-z1(1);
z2 = linspace(zcut,L,n2);
dz2 = z2(2)-z2(1);
tstart = 0;
tend = 1;
dt = 1e-5;
nt = round((tend-tstart)/dt);
t = linspace(tstart,tend,nt);
T_at_zcut = 90;
T_at_0 = 200;
T_at_L = 340;
T = zeros(n1+n2-1,nt);
T(1,:) = T_at_0;
T(end,:) = T_at_L;
T(n1,:) = T_at_zcut;
for i = 1:nt-1
  T(2:n1-1,i+1) = T(2:n1-1,i) + dt/dz1^2*(T(3:n1,i)-2*T(2:n1-1,i)+T(1:n1-2,i));
  T(n1+1:n1+n2-2,i+1) = T(n1+1:n1+n2-2,i) + dt/dz2^2*(T(n1+2:n1+n2-1,i)-...
                          2*T(n1+1:n1+n2-2,i)+T(n1:n1+n2-3,i));
end
plot([z1,z2(2:end)],[T(:,1),T(:,100),T(:,250),T(:,500),T(:,750),T(:,1000),T(:,2500),T(:,5000),T(:,7500),T(:,end)])

[Alfonso]

The following row:

T(n1,:) = T_at_zcut;

Is the one that fix the temperature in the center?

[Torsten]

Yes. The partial differential equation for T is only applied in i=2,...,n1-1,n1+1,...,n1+n2-2. The temperatures in the two boundary points and in the center are kept constant and are excluded from updates over time.

This is equivalent to solving two problems separately:

Problem 1:

dT/dt = d^2T/dz^2   0 <= z <= 0.5
T(z,0) = 0
T(0,t) = 200
T(0.5,t) = 90

Problem 2:

dT/dt = d^2T/dz^2   0.5 <= z <= 1
T(z,0) = 0
T(0.5,t) = 90
T(1,t) = 340

and glue the solutions at the right resp. left boundary points together.

[Alfonso]

Ohh yes, i got this code. What do you think? I think it's working!

function Bioreactor_central

clc
clear all

% Forward Euler method is used for solving the following boundary value problem:

% Problem parameters
alpha1 = 0.005;
mu0 = 5; 
mue = 3.5; 
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2; 
HP = 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 

lambda = 2 ;
rs = 10;

P(1) = 3;
N(1) = 9;
C(1) = 30;

Z = 10;

% Check of the problem's parameters
if  any([P(1)-PC  N(1)-NC] <= 0 ) 
    error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end

% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;

% Initialization
z = linspace(0, Z, nz+1);
t = linspace(0, T, nt+1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt+1);
P = zeros(1, nt+1);
C = zeros(1, nt+1);

I0 = @(t) 100*max(sin((t-0.4*3600)/(6.28)), sin(0));
w = zeros(nz+1, nt+1);
w(50,:) = 5;

% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1+n2-1, nt+1);

% Finite-difference method
for j = 1:nt

    f1 = (KP*(P(j)-PC))/(HP+(P(j)-PC));
    f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC));
    pH = (6.35 - log10(C(j))) /2 ;
    f3 = 1/(1+exp(lambda*(pH-PHs3)));
    f4 = 1/(1+exp(lambda*(pH-PHs4))); 
    
    I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
    g = mu0*I_in./(Hl+I_in);
    
    integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j));
    N(j+1) = N(j) + dt*(-1/Z*integral_term*N(j));
    P(j+1) = P(j) + dt*(-1/Z*integral_term*P(j));
    C(j+1) = C(j) + dt*(-1/Z*integral_term*C(j));
    
    w(nz+1,j+1) = w(nz,j+1);
    w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + alpha1 * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
    w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + alpha1 * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
    w(1,j+1) = w(2,j+1);
    
    plot(z,w(:,j+1),'b');
    xlabel('z'); ylabel('w');
    title(['t=',num2str(t(j)),' s']); 
    axis([0 10 0 10]);
    pause(0.0001);

end
end

Do you think that the results can be associated to the curve that i draw? I think so

[Torsten]

The setting

w(nz+1,j+1) = w(nz,j+1);

seems to be out of order because at the position of the code where this line appears, w(nz,j+1) is still set to 0 from the line of initialization

w = zeros(nz+1, nt+1);

So your boundary conditions for w are

dw/dx(0) = 0, w(L) = 0, w(L/2) = 5

Is this the correct setting you intend ?

[Alfonso]

i would like to have the same boundaries in the two border. I would like to have a neumann but not equal to 0. I would like to have the same behavior of the curve near the borders. I dont want 0.

I mean, in the code there is the light and it has to be put on a single side of the bioreactor, like the side at z=0. While on the other side, at z=10 i don't want the light so that with this plot i will see that from the central value, i have a bigger curve near the side with the light and a smaller curve near the side withtout the light.

[Torsten]

If you want dw/dz = 0 at both ends, use

w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + alpha1 * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + alpha1 * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
w(1,j+1) = w(2,j+1);
w(nz+1,j+1) = w(nz,j+1);

instead of

w(nz+1,j+1) = w(nz,j+1);
w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + alpha1 * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + alpha1 * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
w(1,j+1) = w(2,j+1);

And since you decouple the two regions at z = L/2, I'm not sure if your integrations with trapz and cumtrapz still have to be over the complete region [0 L] or over [0 L/2] and [L/2 L] separately.

[Alfonso]

Ok @Torsten, i modified the code as you told me and so i get this code:

function Bioreactor_central

clc
clear all

% Forward Euler method is used for solving the following boundary value problem:

% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5; 
mue = 2; 
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2; 
HP = 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 

lambda = 2 ;
rs = 10;

P(1) = 3;
N(1) = 9;
C(1) = 30;

Z = 10;

% Check of the problem's parameters
if any([P(1) - PC  N(1) - NC] <= 0 ) 
    error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end

% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;

% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);

I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;

% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);

% Finite-difference method
for j = 1:nt

    f1 = (KP * (P(j) - PC)) / (HP + (P(j) - PC));
    f2 = (KN * (N(j) - NC)) / (HN + (N(j) - NC));
    pH = (6.35 - log10(C(j))) / 2 ;
    f3 = 1 / (1 + exp(lambda * (pH - PHs3)));
    f4 = 1 / (1 + exp(lambda * (pH - PHs4))); 
    
    I_in = I0(t(j)) .* exp(-K * z.') .* exp((-rs) .* cumtrapz(z.', w(:, j)));
    g = mu0 * I_in ./(Hl + I_in);
    
    integral_term = f1 * f2 * f3 * trapz(z.', g .* w(:, j));
    N(j + 1) = N(j) + dt * (-1 / Z * integral_term * N(j));
    P(j + 1) = P(j) + dt * (-1 / Z * integral_term * P(j));
    C(j + 1) = C(j) + dt * (-1 / Z * integral_term * C(j));
    
    w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
    w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
    w(1,j+1) = w(2,j+1);
    w(nz+1,j+1) = w(nz,j+1);
    
    plot(z, w(:, j + 1), 'b');
    xlabel('z'); ylabel('w');
    title(['t=', num2str(t(j)), ' s']); 
    axis([0 10 0 10]);
    pause(0.0001);

end
end

As a plot i have for example:

That seems correct because from the center value, i see 2 different curves with 2 different values on the borders. But, the I_in, which is the irradiation, is it put only on 1 side? Because at least i see that the 2 curves are increasing a lot, yes with 2 different values of w, but they are increasing together and i don't want something like this. In this case the plot should have a curve that has a permanent lower values of W and the other curve which is the oppost, i mean this other curve has to increase the w value with a different velocity due to the fact that the algae are irradiated by the light.

[Torsten]

I cannot help you to check whether your equations are adequate to model what you want to model. I can only help you to check whether predefined equations and boundary conditions are correctly implemented in MATLAB code.

[Alfonso]

Ok, so. The equation of the irradiation is:

I_in = I0(t(j)) .* exp(-K * z.') .* exp((-rs) .* cumtrapz(z.', w(:, j)));

Now, i would like to know, if you can help me, where it is put. I mean i want the irradiation only on 1 side of the rectangule boundary, so it is like this? Or the irradiation is everywhere for all the z values?

[Torsten]

I cannot answer this.

At the moment, I_in is defined for all values of z and enters all equations for all values of z. Maybe you should save and plot the profiles to see what is happening.

[Alfonso]

What do you mean that i should save and plot the profiles, which one, i didn't understand

[Torsten]

From what you write it seems to me that I_in as defined does not fit your needs.

So you could save I_in as you save w in a 2d-array and plot it to see if something goes wrong.

[Alfonso]

@Torsten

I have the following updated code:

function Bioreactor_central

clc
clear all

% Forward Euler method is used for solving the following boundary value problem:

% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5; 
mue = 2; 
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2; 
HP = 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 

lambda = 2 ;
rs = 10;

P(1) = 3;
N(1) = 9;
C(1) = 30;

Z = 10;

% Check of the problem's parameters
if any([P(1) - PC  N(1) - NC] <= 0 ) 
    error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end

% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;

% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);

I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;

% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);

% Finite-difference method
for j = 1:nt

    f1 = (KP * (P(j) - PC)) / (HP + (P(j) - PC));
    f2 = (KN * (N(j) - NC)) / (HN + (N(j) - NC));
    pH = (6.35 - log10(C(j))) / 2 ;
    f3 = 1 / (1 + exp(lambda * (pH - PHs3)));
    f4 = 1 / (1 + exp(lambda * (pH - PHs4))); 
    
    I_in = I0(t(j)) .* exp(-K * z.') .* exp((-rs) .* cumtrapz(z.', w(:, j)));
    g = mu0 * I_in ./(Hl + I_in);
    
    integral_term = f1 * f2 * f3 * trapz(z.', g .* w(:, j));
    N(j + 1) = N(j) + dt * (-1 / Z * integral_term * N(j));
    P(j + 1) = P(j) + dt * (-1 / Z * integral_term * P(j));
    C(j + 1) = C(j) + dt * (-1 / Z * integral_term * C(j));
    
    w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
    w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
    w(1,j+1) = w(2,j+1);
    w(nz+1,j+1) = w(nz,j+1);
    
    plot(z, w(:, j + 1), 'b');
    xlabel('z'); ylabel('w');
    title(['t=', num2str(t(j)), ' s']); 
    axis([0 10 0 10]);
    pause(0.0001);

end
end

I understood that I_in is fixed for every value of the z axis, but instead i would like to have it only for like z=10 for every J value. I've tried to modify it and i got:

function Bioreactor_central

clc
clear all

% Forward Euler method is used for solving the following boundary value problem:

% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5; 
mue = 2; 
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2; 
HP = 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 

lambda = 2 ;
rs = 10;

P(1) = 3;
N(1) = 9;
C(1) = 30;

Z = 10;

% Check of the problem's parameters
if any([P(1) - PC  N(1) - NC] <= 0 ) 
    error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end

% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;

% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);

I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;

% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);

I_in = zeros(nz + 1, nt + 1);
    for j = 1:nt
        I_in(:, j) = I0(t(j)) * exp(-K * z.') .* exp((-rs) * cumtrapz(z, w(:, j).')');
    end

% Finite-difference method
for j = 1:nt

    f1 = (KP * (P(j) - PC)) / (HP + (P(j) - PC));
    f2 = (KN * (N(j) - NC)) / (HN + (N(j) - NC));
    pH = (6.35 - log10(C(j))) / 2 ;
    f3 = 1 / (1 + exp(lambda * (pH - PHs3)));
    f4 = 1 / (1 + exp(lambda * (pH - PHs4))); 
    
    I_in_10 = I_in(10, j);
    g = mu0 * I_in_10 / (Hl + I_in_10);
    
    integral_term = f1 * f2 * f3 * trapz(z, g .* w(:, j));
    N(j + 1) = N(j) + dt * (-1 / Z * integral_term * N(j));
    P(j + 1) = P(j) + dt * (-1 / Z * integral_term * P(j));
    C(j + 1) = C(j) + dt * (-1 / Z * integral_term * C(j));
    
    w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g(2:50-1) + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
    w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g(50+1:nz) + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
    w(1,j+1) = w(2,j+1);
    w(nz+1,j+1) = w(nz,j+1);
    
    plot(z, w(:, j + 1), 'b');
    xlabel('z'); ylabel('w');
    title(['t=', num2str(t(j)), ' s']); 
    axis([0 10 0 10]);
    pause(0.0001);

end
end

The problem is that the code works but doesn't give me anything as result. And i dont know why D:

Can you help me?

[Torsten]

Bioreactor_central()

function Bioreactor_central
% Forward Euler method is used for solving the following boundary value problem:
% Problem parameters
Dm = 5 * 10^(-4);
mu0 = 5; 
mue = 2; 
Hl = 5;
T = 5000;
nz = 100; %nt = 1000;
KP = 2; 
HP = 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 
lambda = 2 ;
rs = 10;
P(1) = 3;
N(1) = 9;
C(1) = 30;
Z = 10;
% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;
% Check of the problem's parameters
if any([P(1) - PC  N(1) - NC] <= 0 ) 
    error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end
% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);
I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;
I_in = zeros(nz + 1,nt);
% Finite-difference method
for j = 1:nt
    f1 = (KP * (P(j) - PC)) / (HP + (P(j) - PC));
    f2 = (KN * (N(j) - NC)) / (HN + (N(j) - NC));
    pH = (6.35 - log10(C(j))) / 2 ;
    f3 = 1 / (1 + exp(lambda * (pH - PHs3)));
    f4 = 1 / (1 + exp(lambda * (pH - PHs4))); 
    I_in(:,j) = I0(t(j)) .* exp(-K * z.') .* exp((-rs) .* cumtrapz(z.', w(:, j)));   
    I_in_10 = I_in(Z,j);
    g_in_10 = mu0 * I_in_10 / (Hl + I_in_10);
    
    integral_term = f1 * f2 * f3 * trapz(z, g_in_10 .* w(:, j));
    N(j + 1) = N(j) + dt * (-1 / Z * integral_term * N(j));
    P(j + 1) = P(j) + dt * (-1 / Z * integral_term * P(j));
    C(j + 1) = C(j) + dt * (-1 / Z * integral_term * C(j));
    
    w(2:50-1,j+1) = w(2:50-1,j) + dt * (f1*f2*f3*w(2:50-1,j).*g_in_10 + Dm * (w(3:50,j)-2*w(2:50-1,j)+w(1:50-2,j))/dz^2);
    w(50+1:nz,j+1) = w(50+1:nz,j) + dt * (f1*f2*f3*w(50+1:nz,j).*g_in_10 + Dm * (w(50+2:nz+1,j)-2*w(50+1:nz,j)+w(50:nz-1,j))/dz^2);
    w(1,j+1) = w(2,j+1);
    w(nz+1,j+1) = w(nz,j+1);
    
end
plot(z,w(:,end))
end

[Alfonso]

@Torsten

Ok, now the code is working again but i'm not having what i would like to have... :C

As you can see, i would like to have the red curve, not the blue which comes from the model. As i understood from the code that you sent the I_in is fixed only on z=10 and it is right but it doesnt happen what it should :/.

The algae should grow faster near the light side then the shadow side.

[Alfonso]

@Torsten

I've tried to split everything and putting I_in equal to 0 from z=0 to z=5, so that i have the following code:

function Bioreactor_central

clc
clear all

% Forward Euler method is used for solving the following boundary value problem:

% Problem parameters
Dm = 5;% * 10^(-4);
mu0 = 5; 
mue = 2; 
Hl = 5;
T = 5000;
nz = 100; nt = 1000;
KP = 2; 
HP = 7; 
PC = 1;
K = 2;
KN = 3;
HN = 14.5 * 10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 

lambda = 2 ;
rs = 10;

P(1) = 3;
N(1) = 9;
C(1) = 30;

P_1(1) = 3;
N_1(1) = 9;
C_1(1) = 30;


Z = 10;

% Check of the problem's parameters
if any([P(1) - PC  N(1) - NC] <= 0 ) 
    error('Check problem parameters')
end
if any([Dm mue mu0 Hl Z T nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end

% Stability
dz = Z / nz;
st = 1 / 2;
dt = st * dz^2 / Dm;
nt = T / dt;

% Initialization
z = linspace(0, Z, nz + 1);
t = linspace(0, T, nt + 1);
nt = round(nt);
nz = round(nz);
N = zeros(1, nt + 1);
P = zeros(1, nt + 1);
C = zeros(1, nt + 1);

N_1 = zeros(1, nt + 1);
P_1 = zeros(1, nt + 1);
C_1 = zeros(1, nt + 1);

I0 = @(t) 100 * max(sin((t - 0.4 * 3600) / (6.28)), sin(0));
w = zeros(nz + 1, nt + 1);
w(50, :) = 5;

% Initialize T with the correct number of columns
n1 = 50;
n2 = 50;
T = zeros(n1 + n2 - 1, nt + 1);


% Finite-difference method
for j = 1:nt

    f1 = (KP * (P(j) - PC)) / (HP + (P(j) - PC));
    f2 = (KN * (N(j) - NC)) / (HN + (N(j) - NC));
    pH = (6.35 - log10(C(j))) / 2 ;
    f3 = 1 / (1 + exp(lambda * (pH - PHs3)));
    f4 = 1 / (1 + exp(lambda * (pH - PHs4))); 

    f11 = (KP * (P_1(j) - PC)) / (HP + (P_1(j) - PC));
    f22 = (KN * (N_1(j) - NC)) / (HN + (N_1(j) - NC));
    pH1 = (6.35 - log10(C_1(j))) / 2 ;
    f33 = 1 / (1 + exp(lambda * (pH1 - PHs3)));
    f44 = 1 / (1 + exp(lambda * (pH1 - PHs4))); 
    
    I_in = zeros(nz + 1, 1);
    I_in_1_5 = 0;%I0(t(j)) .* exp(-K * z(1:5).') .* exp((-rs) .* cumtrapz(z(1:5).', w(1:5, j)));
    I_in_5_10 = I0(t(j)) .* exp(-K * z(5:10).') .* exp((-rs) .* cumtrapz(z(5:10).', w(5:10, j)));
    g_1_5 = mu0 * I_in_1_5 ./(Hl + I_in_1_5);
    g_5_10 = mu0 * I_in_5_10 ./(Hl + I_in_5_10);

    integral_term_1_5 = f11 * f22 * f33 * trapz(z(1:5), g_1_5 .* w(1:5, j));
    integral_term_5_10 = f1 * f2 * f3 * trapz(z(5:10), g_5_10 .* w(5:10, j));
    
    N_1(j + 1) = N(j) + dt * (-1 / Z * integral_term_1_5 * N(j));
    P_1(j + 1) = P(j) + dt * (-1 / Z * integral_term_1_5 * P(j));
    C_1(j + 1) = C(j) + dt * (-1 / Z * integral_term_1_5 * C(j));
    
    N(j + 1) = N(j) + dt * (-1 / Z * integral_term_5_10 * N(j));
    P(j + 1) = P(j) + dt * (-1 / Z * integral_term_5_10 * P(j));
    C(j + 1) = C(j) + dt * (-1 / Z * integral_term_5_10 * C(j));
    
    w(2:49,j+1) = w(2:49,j) + dt * (f11*f22*f33*w(2:49,j).*g_1_5(2:49) + Dm * (w(3:50,j) - 2*w(2:49,j) + w(1:48,j))/dz^2) - mue*f4*w(2:49,j);
    w(51:nz,j+1) = w(51:nz,j) + dt * (f1*f2*f3*w(51:nz,j).*g_5_10(51:nz) + Dm * (w(52:nz+1,j)-2*w(51:nz,j)+w(50:nz-1,j))/dz^2) - mue*f4*w(51:nz,j);

    w(1,j+1) = w(2,j+1);
    w(nz+1,j+1) = w(nz,j+1);
    
    plot(z, w(:, j + 1), 'b');
    xlabel('z'); ylabel('w');
    title(['t=', num2str(t(j)), ' s']); 
    axis([0 10 0 10]);
    pause(0.0001);

end
%plot(z,w(end));
end

But i have the following error:

Index exceeds the number of array elements. Index must not exceed 1.

Error in Bioreactor_central1 (line 110)
    w(2:49,j+1) = w(2:49,j) + dt * (f11*f22*f33*w(2:49,j).*g_1_5(2:49) + Dm * (w(3:50,j) - 2*w(2:49,j) + w(1:48,j))/dz^2) - mue*f4*w(2:49,j);

Do you know how can i solve?

[Torsten]

You split the domain and decouple all the equations over z=5. So why not making two computations on 0 <= z <= 5, once with I_in = 0 and second with the usual I_in ?

[Alfonso]

@Torsten

Yes, that's what i would like to do but it's hard for me to think about what to modify on matlab to obtain what i want xD

[Torsten]

Use your old code without the splitting of the w-array and the length of the region halfed - that's all.

[Alfonso]

You are perfectly right! I've returned to the old code, the following:

function w = Bioreattore

clc
clear all

% Forward Euler method is used for solving the following boundary value problem:

% Problem parameters
alpha1 =  5*10^(-4); 
mu0 = 5; 
mue = 3.5; 
Hl = 5;
T = 10000;                                                  % [s]
nz = 100; nt = 1000;
KP = 2 ; 
HP= 7; 
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 

lambda = 2 ; 
rs= 10;

P(1)=3;                                                     % Initial chosen concentration of P
N(1)=9;                                                     % Initial chosen concentration of N 
C(1)=30;                                                    % Initial chosen concentration of C

Z = 10;                                                     % Lenght of the bioreactor

%Check of the problem's parameters
if  any([P(1)-PC  N(1)-NC] <=0 ) 
    error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T  nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end

% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;

% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt);                                             % Round nt to an integer
nz = round(nz);                                             % Round nz to an integer
N = zeros(1, nt+1);                                         % Initialize N to zero
P = zeros(1, nt+1);                                         % Initialize P to zero
C = zeros(1, nt+1);                                         % Initialize C to zero

I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w = zeros(nz+1,nt+1);
w(:,1) = 0;                                                 %Final concentration of the algae


% Finite-difference method
for j = 1:nt

    f1= (KP*(P(j)-PC))/(HP+(P(j)-PC));
    f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC));
    pH = (6.35 - log10(C(j))) /2 ;
    f3 = 1/(1+exp(lambda*(pH-PHs3)));
    f4 = 1/(1+exp(lambda*(pH-PHs4))); 
    
    I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
    g = mu0*I_in./(Hl+I_in);
    
    
    integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j));    % Data organization
    N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j));           % Function of N over the time and space
    P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j));           % Function of P over the time and space
    C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j));           % Function of C over the time and space
    
    w(nz+1,j+1) = 5;                                        % Algae concentration at z=0    
    w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + alpha1 * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2)- mue*f4*w(2:nz,j);
    w(1,j+1) = w(2,j+1);                                    % Algae concentration at z=10
   
    
    plot(z,flipud(w(:,j+1)),'b');
    %set(gca,'XDir','reverse');
    xlabel('z'); ylabel('w');
    title(['t=',num2str(t(j))]); 
    axis([0 10 0 10]);
    pause(0.0001);
end
end

But at least while i remove put I_in equal to 0 i have the same plot, as follows:

What it's you advise?

[Torsten]

My advice remains the same:

Plot

I_in=I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j))); 

as you plot w to see if it's equal to or almost equal to 0.

[Alfonso]

@Torsten

Ok, so. I did it, the code is the following:

function w = Bioreattore_Matlab 
clc
clear all
% Forward Euler method is used for solving the following boundary value problem:
% U(0,t) = 0
% U(L,t) = u (end-1,t)
% U(z, 0) = U0(z)=0 
% it is assumed that the concentration of C,N,P are costant
% Problem parameters
alpha1 =  5*10^(-4); 
mu0 = 5; 
mue = 3.5; 
Hl = 5;
T = 10000; % [s]
nz = 100; nt = 1000;
KP = 2 ; 
HP= 7; 
%P = 3;
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
%N = 4;
NC = 2;
KC = 5;
HC = 5 ;
%C = 3 ;
PHs3 = 4 ;
PHs4 = 5; 
Z = 10;
lambda = 2 ; 
rs= 10;

P(1)=3;
N(1)=9;
C(1)=30;
%Check of the problem's parameters
if  any([P(1)-PC  N(1)-NC] <=0 ) 
    error('Check problem parameters')
end
if any([alpha1 mu0 Hl Z T  nz-2 nt-2 KP P PC HC KN N NC KC HC C HN ] <= 0)
    error('Check problem parameters')
end
% Stability
dz=Z/nz ;
%dt = T/nt;
%r = alpha1*dt/((dz)^2);
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;
%while (2*r) > st
 %   nt = nt+1;
%end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt); % Arrotonda nt a un numero intero
N = zeros(1, nt+1); % Inizializzazione di N a zero
P = zeros(1, nt+1); % Inizializzazione di P a zero
C = zeros(1, nt+1); % Inizializzazione di C a zero

%f1= (KP*(P-PC))/(HP+(P-PC));
%f2 = (KN*(N-NC))/(HN+(N-NC));
%pH = (6.35 - log10(C)) /2 ;
%f3 = 1/(1+exp(lambda*(pH-PHs3)));
%f4 = 1/(1+exp(lambda*(pH-PHs4))); 

%I0 = 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
nz = round(nz);
nt = round(nt);
w = zeros(nz+1,nt+1);
w(:,1) = 0; %Final concentration of the algae [0]
% Finite-difference method
%for j=2:nt+1
%   u(1) = 0; % Condizione al contorno di Dirichlet
%   g = (mu0.*I(j))/(Hl+I(j)); decay = g(j)*f1*f2*f3-mue*f4;
%  u(2:end-1)=(1-2*r-decay*dt)*u(2:end-1)+r*(u(1:end-2)+u(3:end));
%   u(end) = u(end); % Condizione al contorno di Neumann
for j = 1:nt

    f1= (KP*(P(j)-PC))/(HP+(P(j)-PC));
    f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC));
    pH = (6.35 - log10(C(j))) /2 ;
    f3 = 1/(1+exp(lambda*(pH-PHs3)));
    f4 = 1/(1+exp(lambda*(pH-PHs4))); 
    
    I_in = I0(t(j)).*exp(-K*z.').*exp((-rs));%.*cumtrapz(z.',w(:,j)));
    g = mu0*I_in./(Hl+I_in);
    integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j));
    
    N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j))
    P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j))
    C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j))
    
    
    %I_in = I0(t(j))*exp(-K*z.').*exp(-rs)*cumtrapz(z.',w(:,j));
    
    
        
    %w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz,j)+r*(w(1:nz-1,j)+w(3:nz+1,j));
    w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + alpha1 * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2)
    w(1,j+1) = w(2,j+1);
    w(nz+1,j+1) = 5;%w(nz,j+1); %Concentration at z=0
    %w(2:nz,j+1) = w(2:nz,j) + dt*...
    
    %w(2:nz,j+1) = integral_term(j) +  
    
    plot(z,flipud(I_in),'b');
    %set(gca,'XDir','reverse');
    xlabel('z'); ylabel('w');
    title(['t=',num2str(t(j))]); 
    %axis([0 10 0 10]);
    pause(0.0001);
end
%plot(z,flipud(w(:,end)))
end

But the plots are not good as i expected, look at that:

They are the only plot that i see for all the time t is like 1 second the first plot and the next second the other. What do you think? What's you advise, they are just the I_in plot, dont read on yaxes the label.

[Alfonso]

@Torsten

Regarding the old code:

function w = Bioreattore

clc
clear all

% Forward Euler method is used for solving the following boundary value problem:

% Problem parameters
alpha1 =  5*10^(-4); 
mu0 = 0; 
mue = 99999999;%3.5; 
Hl = 5;
T = 15000;                                                  % [s]
nz = 100; nt = 1000;
KP = 2 ; 
HP= 7; 
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
NC = 2;
KC = 5;
HC = 5 ;
PHs3 = 4 ;
PHs4 = 5; 

lambda = 2 ; 
rs= 10;

P(1)=3;                                                     % Initial chosen concentration of P
N(1)=9;                                                     % Initial chosen concentration of N 
C(1)=30;                                                    % Initial chosen concentration of C

Z = 10;                                                     % Lenght of the bioreactor

%Check of the problem's parameters
%if  any([P(1)-PC  N(1)-NC] <=0 ) 
%    error('Check problem parameters')
%end
%if any([alpha1 mu0 Hl Z T  nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
%    error('Check problem parameters')
%end

% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/alpha1;
nt = T/dt;

% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt);                                             % Round nt to an integer
nz = round(nz);                                             % Round nz to an integer
N = zeros(1, nt+1);                                         % Initialize N to zero
P = zeros(1, nt+1);                                         % Initialize P to zero
C = zeros(1, nt+1);                                         % Initialize C to zero

I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w = zeros(nz+1,nt+1);
w(:,1) = 0;                                                 %Final concentration of the algae


% Finite-difference method
for j = 1:nt

    f1= (KP*(P(j)-PC))/(HP+(P(j)-PC));
    f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC));
    pH = (6.35 - log10(C(j))) /2 ;
    f3 = 1/(1+exp(lambda*(pH-PHs3)));
    f4 = 1/(1+exp(lambda*(pH-PHs4))); 
    
    I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
    g = mu0*I_in./(Hl+I_in);
    
    
    integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j));     % Data organization
    N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j));           % Function of N over the time and space
    P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j));           % Function of P over the time and space
    C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j));           % Function of C over the time and space
    
    w(nz+1,j+1) = 5;                                         % Algae concentration at z=0    
    w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + alpha1 * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2)- mue*f4*w(2:nz,j);
    w(1,j+1) = w(2,j+1);                                     % Algae concentration at z=10
   
    
    plot(z,flipud(w(:,j+1)),'b');
    %set(gca,'XDir','reverse');
    xlabel('z'); ylabel('w');
    title(['t=',num2str(t(j))]); 
    axis([0 10 0 10]);
    pause(0.0001);
end
end

I see that if i change every parameter it doesn't change anything in the curve, itìs not normal, it seems that the w equation doesnt read the initial parameter, it gives me everytime the same curve also if i put everything 0. I just need a code that reads those changes at least, can you check?

[Torsten]

The parameters are all correctly used, but they are most probably far off. Did you take care about the units? In the article, some of them have time unit "day" and others have time unit "second". Did you convert them correctly to a common unit in your code ?

And the setting

w(nz+1,j+1) = 5;

is at the wrong position.

It should be set before the for-loop as

w(nz+1,:) = 5

[Alfonso]

Ok, man i will check. At least i didn't check the dimension before D:

[Alfonso]

@Torsten also with the right dimension the code doesn't work very well. Can i ask you a question? If i have this code:

function w = Bioreattore

clc
clear all

% Forward Euler method is used for solving the following boundary value problem:

% Problem parameters
Dm =  5*10^(-4);                                            %m^2/s              
mu0 = 8000%2.5;                                                  % L/day
mue = 3.5;                                                  % 1/day

Hl = 5;     %half-saturation intensity W/(m^2*day)
T = 10000;                                                  % [s]
nz = 100; nt = 1000;
KP = 2 ; 
HP= 7;             %half-saturation concentrations of Phosphates
PC = 1;     %critical concentration of the phosphates dopo la quale la crescita di w va a 0.
K = 2;
KN = 3 ;
HN = 0.02;   %14.5*10^(-6);  %half-saturation concentrations of nitrates [mgN/L]
NC = 2;     %critical concentration of the nitrates dopo la quale la crescita di w va a 0.
KC = 5;
HC = 5 ;
PHs3 = 9 ; %describes the "switching" value of pH at which the growth increases if all other factors are kept unchanged.
 

lambda = 2 ;                %sharpness of the profile of the f3 function
rs= 10;                                                     %L*m/d

P(1)=50;                                                     % Initial chosen concentration of P [mg/L]
N(1)=330;                                                     % Initial chosen concentration of N [mg/L]
C(1)=2100;                                                    % Initial chosen concentration of C/ Aumentando C, si abbassa il pH(+acido) e aumenta la morte algale [mg/L]

Z = 10;                                                     % Lenght of the bioreactor

%Check of the problem's parameters
if  any([P(1)-PC  N(1)-NC] <=0 ) 
    error('Check problem parameters')
end
if any([Dm mu0 Hl Z T  nz-2 nt-2 KP P(1) PC HC KN N(1) NC KC HC C(1) HN ] <= 0)
    error('Check problem parameters')
end

% Stability
dz = Z/nz;
st = 1/2;
dt = st*dz^2/Dm;
nt = T/dt;

% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
nt = round(nt);                                             % Round nt to an integer
nz = round(nz);                                             % Round nz to an integer
N = zeros(1, nt+1);                                         % Initialize N to zero
P = zeros(1, nt+1);                                         % Initialize P to zero
C = zeros(1, nt+1);                                         % Initialize C to zero

I0 = @(t) 100*max(sin ((t-0.4*3600)/(6.28)) , sin (0));
w = zeros(nz+1,nt+1);
w(:,1) = 0;                                                 %Final concentration of the algae


% Finite-difference method
for j = 1:nt

    f1= (KP*(P(j)-PC))/(HP+(P(j)-PC))
    f2 = (KN*(N(j)-NC))/(HN+(N(j)-NC))
    pH = (6.35 - log10(C(j))) /2   %relation between the pH and the concentration of carbon dioxide
    f3 = 1/(1+exp(lambda*(pH-PHs3)))   %The growth rate dependence function. Si riduce se il tutto diventa più acido (pH scende)..
    
    Ha = mue*f3;       %death rate of the algae. [gCOD/L]
    
    I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));     % W/m^2*s
    g = mu0*I_in./(Hl+I_in);
    
    
    integral_term = f1 * f2 * f3 * trapz(z.',g.*w(:,j));     % Data organization
    N(j+1) = N(j) + dt*( -1/Z*integral_term*N(j));           % Function of N over the time and space
    P(j+1) = P(j) + dt*( -1/Z*integral_term*P(j));           % Function of P over the time and space
    C(j+1) = C(j) + dt*( -1/Z*integral_term*C(j));           % Function of C over the time and space
    
    w(nz+1,j+1) = 5;                                         % Algae concentration at z=0    [gSST/L]
    w(2:nz,j+1) = w(2:nz,j) + dt*(f1*f2*f3*w(2:nz,j).*g(2:nz) + Dm * (w(3:nz+1,j)-2*w(2:nz,j)+w(1:nz-1,j))/dz^2)- Ha*w(2:nz,j);
    w(1,j+1) = w(2,j+1);                                     % Algae concentration at z=10
   
    
    plot(z,flipud(w(:,j+1)),'b');
    %set(gca,'XDir','reverse');
    xlabel('z'); ylabel('w');
    title(['t=',num2str(t(j))]); 
    axis([0 10 0 10]);
    pause(0.0001);
end
end

How can i modify the plot to obtain the following graphs: