Unable to reproduce a working code for threepointbvp

Hi everyone,

I am struggling to reproduce a working threepointbvp code. Please see attached or below. What am I doing wrong?

The working code is:

function threebvp(solver)
%THREEBVP  Three-point boundary value problem
%   This example of multipoint BVPs comes from the study of a physiological
%   flow problem described in C. Lin and L. Segel, Mathematics Applied to
%   Deterministic Problems in the Natural Sciences, SIAM, Philadelphia, PA,
%   1988. For x in [0, lambda], the equations for the problem are
%
%       v' = (C - 1)/n
%       C' = (vC - min(x,1))/eta
%
%   for known parameters n, kappa, and lambda > 1. eta = lambda^2/(n*kappa^2).
%   The term min(x,1) in the equation for C'(x) is not smooth at x = 1.
%   Lin and Segel describe this BVP as two problems, one set on [0, 1]
%   and the other on [1, lambda], connected by the requirement that v(x)
%   and C(x) be continuous at x = 1. The solution is to satisfy the boundary
%   conditions v(0) = 0 and C(lambda) = 1.  BVP4C solves this problem as
%   a three-point BVP, imposing internal conditions at the interface point
%   x = 1.
%
%   This example solves the problem for n = 5e-2, lambda = 2, and a range
%   kappa = 2,3,4,5. The solution for one value of kappa is used as guess
%   for the next.
%
%   By default, this example uses the BVP4C solver. Use syntax
%   THREEBVP('bvp5c') to solve this problem with the BVP5C solver.
%
%   See also BVP4C, BVP5C, BVPSET, BVPGET, BVPINIT, DEVAL, FUNCTION_HANDLE.

%   Jacek Kierzenka and Lawrence F. Shampine


if nargin < 1
   solver = 'bvp4c';
end
bvpsolver = fcnchk(solver);

% Problem parameter, shared with nested functions
n = 5e-2;

% Initial mesh - duplicate the interface point x = 1.
lambda = 2;
xinit = [0, 0.25, 0.5, 0.75, 1, 1, 1.25, 1.5, 1.75, 2];

% Constant initial guess for the solution
yinit = [1; 1];

% The initial profile
sol = bvpinit(xinit,yinit);

% For each kappa, the quantity of interest is the emergent osmolarity.
% This example compares the value computed by BVP4C, Os = 1/v(lambda),
% with Lin and Segel's asymptotic approximation.
fprintf(' kappa    computed Os    approximate Os \n')
for kappa = 2:5
   eta = lambda^2/(n*kappa^2);  % use in nested functions
   sol = bvpsolver(@f,@bc,sol);

     K2 = lambda*sinh(kappa/lambda)/(kappa*cosh(kappa));
     approx = 1/(1 - K2);
     computed = 1/sol.y(1,end);
     fprintf('  %2i    %10.3f     %10.3f \n',kappa,computed,approx);
  end

figure
plot(sol.x,sol.y(1,:),sol.x,sol.y(2,:),'--')
legend('v(x)','C(x)')
title(['A three-point BVP solved with ',upper(func2str(bvpsolver))])
xlabel(['\lambda = ',num2str(lambda),', \kappa = ',num2str(kappa),'.'])
ylabel('v and C')

% -----------------------------------------------------------------------
% Nested functions -- n and eta are provided by the outer function.
%

     function dydx = f(x,y,region)
        % Derivative function -- share n and eta with the outer function.
        dydx = zeros(2,1);
        dydx(1) = (y(2) - 1)/n;
        % The definition of C'(x) depends on the region.
        switch region
           case 1    % x in [0 1]
              dydx(2) = (y(1)*y(2) - x)/eta;
           case 2    % x in [1 lambda]
              dydx(2) = (y(1)*y(2) - 1)/eta;
           otherwise
              error('MATLAB:threebvp:BadRegionIndex','Incorrect region index: %d',region);
        end
     end
  % -----------------------------------------------------------------------

     function res = bc(YL,YR)
        % Boundary (and internal) conditions
        res = [ YL(1,1)             % v(0) = 0
           YR(1,1) - YL(1,2)   % continuity of v(x) at x = 1
           YR(2,1) - YL(2,2)   % continuity of C(x) at x = 1
           YR(2,end) - 1    ]; % C(lambda) = 1
     end
  % -----------------------------------------------------------------------

end  % threebvp

where as my code is:

function threepointbvp(solver)

global A B C D E F G H I J K L M N O P Q R S T U V

%PARAMETERS
A = sym(9.4e-4);     % ks = solid-liquid mass transfer (m/s)
B = sym(3.69e2);     % Ap = solid-liquid surface area (m2/m3)
C = sym(2.89e-9);    % DSO2 = diffusivity of SO2 (m2/s)
D = sym(1.6e-9);     % DCO32 = diffusivity of CO32- (m2/s)
E = sym(4.92e-5);    % CBs = concentration of Ca2+ on limestone surface (mol/m3) 
F = sym(1.6e-8);     % DH = diffusivity of H+ (m2/s)
G = sym(6.24);       % KSO2 = SO2 dissociation equilibrium (mol/m3)
H = sym(2.35e-9);    % DHSO3 = diffusivity of HSO3 (m2/s)
I = sym(5.68e-5);    % KHSO3 = HSO3 dissociation equilibrium constant (mol/m3)
J = sym(1.69e-9);    % DSO32 = diffusivity of SO32 (m2/s)
K = sym(149);        % HSO2 = Henry's constant (m3Pa/mol)
L = sym(47.28325);   % PSO2 = partial pressure of SO2 (Pa) 
M = sym(4.14e-6);    % kga = gas-side mass transfer product (1/s)
N = sym(8.314);      % R = universal gas constant (J/mol.K)
O = sym(323.15);     % T = reactor tempeature (K)
P = sym(2.09);       % CH = concentration of H+ at the gas-liquid interface (mol/m3)
Q = sym(9.333E-06);  % CH = concentration of H+ at the in the bulk slurry (mol/m3)
R = sym(0.679);      % S(IV) = total sulfur concentration in the bulk slurry (mol/m3)
S = sym(4.88e-6);    % C(IV) = total sulfur concentration in the bulk slurry (mol/m3)
T = sym(1.7e-3);     % KCO2 = CO2 dissociation equilibrium constant (mol/m3)
U = 6.55e-8;         % KHCO3 = HCO3 dissociation equilibrium constant (mol/m3)
V = 2.08e-9;         % DHCO3 = Diffusivity of HCO3 (m2/s)

% VARIABLES

% y(1) = CSO2
% y(2) = CHSO3
% y(3) = CSO32
% y(4) = CHCO32
% y(5) = CCO32

% EQUATIONS

% Region I

% dy(1)/dx = y(6);
% d^2y(1)/dx^2 = (2*A*B / C)*(1 + (C*y(1) / 2*D*E) + (F*G*(y(2)/y(1)) / D*E)*E;
% dy(2)/dx = y(7);
% d^2y(2)/dx^2 = -(2*A*B / C)*(1 + (C*y(1) / 2*D*E) + (F*G*(y(2)/y(1)) /D*E)*E;

%Region II

% dy(3)/dx = y(8)
% dy(4)/dx = y(4);
% d^2y(3)/dx^2 = -(A*B / D)*(1 + (F*I*(y(3) / y(2)))/(D*E))*E;
% dy(5)/dx = y(5)

if nargin < 1
   solver = 'bvp4c';
end

bvpsolver = fcnchk(solver);

xinit = [0.05, 0.10, 0.15, 0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.5,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.00,];
yinit = [4.03e-08, 3.39e-02, 2.60e-01, 5.45e-01, 4.82e-03, 0, 0, 0];

solinit = bvpinit(xinit,yinit);
sol = bvpsolver(@materialbalanceodes,@materialbalancebc,solinit);

x = linspace(0.05,1,0.05);
y = deval(sol,x);
figure
plot(x,y(1,:), x,y(2,:), x,y(3,:), x,y(4,:), x,y(5,:), x,y(6,:), x,y(7,:), x,y(8,:))
legend('SO_2','HSO_{3}^{-}','SO_{3}^{2-}','HCO_{3}^{-}','CO_{3}^{2-}')

for x = 0.05:0.5 
    region = 'case 1';    
    for x = 0.55:1
    region = 'case 2';

      end

end

% -----------------------------------------------------------------------

function dydx = materialbalanceodes(x,y, region)
dydx = zeros(8,1);

        switch region
           case 1    % x in [0 0.5]
        dydx(1) = y(6);
        dydx(2) = (2 * A * B / C) * (1 + ((C * y(1)) / (2 * D * E)) + ((F * G )* (y(2)/y(1)) / (D*E)))*E;
        dydx(3) = y(7);
        dydx(4) = -((2 * A * B) / (H)) * (1 + ((C * y(1)) / (2 * D * E)) + ((F * G) * (y(2) / y(1)) / (D * E))) * E;

            case 2    % x in [0.5 1]
        dydx(5) = y(8);
        dydx(6) = y(4);
        dydx(7) = -((A * B) / D) * (1 + (F * I * (y(3) / y(2)))/(D * E)) * E;
        dydx(8) = y(5);  
           otherwise
        error('MATLAB:threebvp:BadRegionIndex','Incorrect region index: %d',region);

end

% -----------------------------------------------------------------------------------------

     function res = materialbalancebc(ya,yb)

        res = [ (ya(1,0.05) - (K * L))                                                    % x = 0.05                          
            (ya(2,0.05) - ((I * K * L) / P)) 
            (ya(6,0.05) +  M * ((L /( N * O))- (K * L)))           
                  ya(1,0.5)                                                               % x = 0.5
                  ya(2,0.5) - ((R*G*Q)/(Q.^2 + G * Q + G * I))
                  ya(3,0.5)
                  ya(4,0.5)
                  - C * ya(6,0.5) - V * ya(4,0.5) - J * yb(8,0.5)
                  - H * yb(7,0.5) + H * yb(7,0.5) - V * ya(4,0.5) - 2 * J * ya(8,0.5)        
                                          (yb(2,1) - ((R * G * Q)/(Q.^2 + G * Q + G * I)))    % x = 1       
                                          (yb(3,1) - ((R * G * I)/(Q.^2 + G * Q + G * I)))
                                          yb(4,1) - ((S * T * Q)/(Q.^2 + T * Q + T * U))]

     end
  % ----------------------------------------------------------------------------------------------

end
end

Where am I going wrong?

Copyright 1984-2014 The MathWorks, Inc.

[Rik]

Today I formatted your code, next time, use the {}Code button. See here for a GIF showing how to do it.

[Rik]

If it is working, why are you changing it? Also, why are you using global variables? They should be avoided if possible, or else you should have long names that will avoid any name collision (which would cause bugs that are very difficult to solve).

[Dursman Mchabe]

Hi Rik, thanks a lot for the response. The reason I'm changing it is that I must use my equations and boundary conditions.