Im working with algae growth linked to the light. But i have some problems on my code

My code is the following:

function w = Bioreattore_POTENTE_(Mat)  
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 = 1000; % [s]
nz = 100; nt = 100;
KP = 2 ; 
HP= 7; 
P = 3 ;
PC = 1;
K = 2;
KN = 3 ;
HN = 14.5*10^(-6);
N = 3;
NC = 2;
KC = 5;
HC = 5 ;
C = 4 ;
PHs3 = 4 ;
PHs4 = 5; 
Z = 10;
lambda = 2 ; 
rs= 10;
%Check of the problem's parameters
if  any([P-PC  N-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;
while (2*r)  > st
nt = nt+1;
end
% Initialization
z = linspace(0,Z,nz+1); t = linspace(0,T,nt+1);
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));
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
    I_in = I0(t(j)).*exp(-K*z.').*exp((-rs).*cumtrapz(z.',w(:,j)));
    %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));
    %w(1,j+1) = w(2,j+1);
    %w(2:nz,j+1)= (1-2*r-integral_term*dt)*w(2:nz-1,j)+r*(w(2:nz-3,j)+w(2:nz-2,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(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,w(:,j+1),'r*:');
    %set(gca,'XDir','reverse');
    xlabel('z'); ylabel('w');
    title(['t=',num2str(t(j))]); 
    axis([0 10 0 10]);
    pause(0.0001);
end
%plot(z,w(:,end))
%end

Now, first of all my code works good until the time T=1000, after this value the code stops working and i don't know why. Second, the script starts at z=10 with a concentration value that i've chosen thanks to the following row in the for cycle:

w(1,j+1) = w(2,j+1); %Concentration at Z=0
w(nz+1,j+1) = 5; %Concentration at Z=10

But i want that the script starts at z=0 with a concentration value that i've chosen. The problem is that when i write the following row in the flow cycle:

w(1,j+1) = 5; %w(2,j+1); Concentration at Z=0
w(nz+1,j+1) = w(nz,j+1); %5; Concentration at Z=10

The script work but the plot not. I have the following, respectively, plots:

z=10; w=5

z=0; w=5

Do you have some suggestions to solve my ploblems? Thank you very much

As attachment you can find the file that im following.

[Image Analyst]

What are you passing in for mat?

See this link to track down the error:

Debugging in MATLAB | Doug's MATLAB Video Tutorials

Answer 1

[Ganesh]

Hi @Alfonso,

The issue you are facing is due to the following lines:

% At T=1001, the loop conditions are met and results into an infinite loop
while (2*r)  > st
    nt = nt+1;
end

The while statement has an update condition that doesn't seem to modify the loop condition. I understand you may need to calculate new values of `r` based on the updated value of nt, for which you may need to use a make-do "do-while loop".

Kindly refer to the Answer to know more on implementing the same: https://in.mathworks.com/matlabcentral/answers/115403-do-while-loop-in-matlab

While this serves as the answer to your question, it would be a better option for you to learn the debugging techniques as mentioned by @Image Analyst.

Kindly refer to the following documentation to know more on debugging MATLAB Code Files:

https://in.mathworks.com/help/matlab/matlab_prog/debugging-process-and-features.html

Hope this helps!