Error in running the code for bvp4c with second order equation, how can I solve this issue?

I will try to explain briefly what happened, since i can't identify the location of my error.

I typed in my 2nd order equation in the "bvpode5" script as follows;

function dydx = bvpode5(x,y)

global Pe Theta K A ;
%
dydx = [y(2) 
  Pe(y(2) -Theta*((y(1)*(A-y(1)))/(K*(A-y(1))+y(1))))];
end

while the original equation is like this 1/Pe (d^2 α)/(dx^2 ) - dα/dx + θ[( α)/(K(A-α)+α)][A-α] = 0

my y(1) is α , while y(2) is dα/dx

As for my boundary condition's script is as follows:

function res = bvpbc5(ya,yb)
global Pe
res = [ya(1)-1-((1/Pe)*ya(2))
    yb(2)
    ];

end 

ya(1) represents α = 1 + (1/Pe)(dα/dx) at x = 0

ya(2) represents dα/dx = 0 at x = 1

Finally, the main script is as follows:

global Pe   Theta K  A;
% 

 % Kc = 0.482;  %   gCOD/gVSS
 Theta = 40;  %  day/day
 % Y = 0.212;   %   gVSS/gCOD
 Pe = 1000; 
 K = 0.123; 
 % Xo = 0.1;  % g/L
 % So = 0.1;  % g/L 
 A= 5.72;   % A = ((Xo/(Y*So))+1)
 xlow = 0;
 xhigh = 1; 
 y=[];
 solinit = bvpinit(linspace(xlow,xhigh,100),[1 0]);
 sol = bvp4c(@bvpode5,@bvpbc5,solinit);
 plot(sol.x,sol.y(1,:))
 xlabel('Reactor Length Z/L')
 ylabel('Substrate Fraction')
 axis ([0 1 0.2 1]);
%

fh=figure(1);
set(fh,'color','white')
title('Substrate Profile Across the Reactor')
disp (sol.y(1,end));

I tried to run the code from this script. However, it shows me this error, how can I solve this issue?;

MainpADM
Array indices must be positive integers or logical values.

Error in bvpode5 (line 6)
  Pe(y(2)-Theta*((y(1)*(A-y(1)))/(K*(A-y(1))+y(1))))];

Error in bvparguments (line 105)
    testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 130)
    bvparguments(solver_name,ode,bc,solinit,options,varargin);

Error in MainpADM (line 16)
 sol = bvp4c(@bvpode5,@bvpbc5,solinit);
 

Answer 1

[Star Strider]

There is a missing operator (most likely multiplication):

Pe(y(2) -Theta*((y(1)*(A-y(1)))/(K*(A-y(1))+y(1))))];
  ↑ ← HERE

This will probably work:

  Pe*(y(2) -Theta*((y(1)*(A-y(1)))/(K*(A-y(1))+y(1))))];

That should at least solve the problem you posted. I did not see any other problems with your code, however I did not run it.

It is always best to not use global variables. Pass the other parameters as additional parameters instead. See Passing Extra Parameters for an extended discussion.

[Hussein Al-Tamreh]

Thank you very much for pointing this out, it is funny that simple things slip through us sometimes.

However, when I plot the graph it shoud be decreasing, but in this code it increases, is the problem with my main equation for defining 2nd order?

according to my reference paper, it should start from the top then decreasing to a value of (1,1) or in between.

Reference graph:

Generated graph:

[Star Strider]

I am not certain what the problem is, however ‘sol.y(2,:)’ appears to be closer to the solution you want, instead of ‘sol.y(1,:)’.

Other than that, check the constants (parameters) and be certain the expreession is coded correctly, since it otherwise appears to do what you want.

The (slightly revised) code —

function dydx = bvpode5(x,y,Pe,Theta,K,A)
%
dydx = [y(2) 
        Pe*(y(2) -Theta*((y(1)*(A-y(1)))/(K*(A-y(1))+y(1))))];
end

function res = bvpbc5(ya,yb,Pe)
res = [ya(1)-1-((1/Pe)*ya(2))
       yb(2)];
end 

guess = @(x) [exp(x); -exp(-x)];

 % Kc = 0.482;  %   gCOD/gVSS
 Theta = 40;  %  day/day
 % Y = 0.212;   %   gVSS/gCOD
 Pe = 1000; 
 K = 0.123; 
 % Xo = 0.1;  % g/L
 % So = 0.1;  % g/L 
 A= 5.72;   % A = ((Xo/(Y*So))+1)
 xlow = 0;
 xhigh = 1; 
%  y=[];
 solinit = bvpinit(linspace(xlow,xhigh,150),guess);
 sol = bvp4c(@(x,y)bvpode5(x,y,Pe,Theta,K,A), @(ya,yb)bvpbc5(ya,yb,Pe),solinit);
  
 figure
 plot(sol.x,sol.y(2,:))
 xlabel('Reactor Length Z/L')
 ylabel('Substrate Fraction')
%  axis ([0 1 0.2 1]);
axis([0  1    -1  10])
%
fh=figure(1);
% set(fh,'color','white')
title('Substrate Profile Across the Reactor')
disp (sol.y(1,end));

.

[Hussein Al-Tamreh]

I will explain what is the objective here, I have to find the behavior of my Axial Dispersion Model (ADM) by defining the dimensionless equation as shown below:

and using bvp4c, I want to plot the Substrate fraction (α) VS Reactor length (x) or (Z/L) {they are all dimensionless}

It looks wonderful, but the Substrate fraction is too high, and the curve is too steep. I think it should be smoothly curved but should not reach (or approach zero) since this reactor is not an ideal case (like Plug Flow Reactor). I have to compare it to the other model to see if it is right or wrong.

I changed the constants' values to see what changes could occur, it seems that Substrate fraction will always approach zero. is it becuase of this code;

guess = @(x) [exp(x); -exp(-x)];

[Star Strider]

Change the ‘guess’ function to something that works for you. Everything else I tried failed, including the original [0 1] you tried.

Meanwhile, I actually did solve your original problem.

I will delete my Answer in a few hours.

[Hussein Al-Tamreh]

Your answers are pretty useful for me, thank you very much in helping me.

[Star Strider]

As always, my pleasure!