This example shows how to solve a multipoint boundary value problem, where the solution of interest satisfies conditions inside the interval of integration.

For $\mathit{x}$ in $\left[0,\lambda \right]$, consider the equations

${\mathit{v}}^{\prime}=\frac{\mathit{C}-1}{\mathit{n}}$,

${\mathit{C}}^{\prime}=\frac{\mathrm{vC}-\mathrm{min}\left(\mathit{x},1\right)}{\eta}$.

The known parameters of the problem are $\mathit{n}$, $\kappa $, $\lambda >1$, and $\eta =\frac{{\lambda}^{2}}{\mathit{n}\cdot {\kappa}^{2}}$.

The term $\mathrm{min}\left(\mathit{x},1\right)$ in the equation for ${\mathit{C}}^{\prime}\left(\mathit{x}\right)$ is not smooth at $\mathit{x}=1$, so the problem cannot be solved directly. Instead, you can break the problem into two: one set in the interval $\left[0,1\right]$, and the other set in the interval $\left[1,\lambda \right]$. The connection between the two regions is that the solutions must be continuous at $\mathit{x}=1$. The solution must also satisfy the boundary conditions

$\mathit{v}\left(0\right)=0$,

$\mathit{C}\left(\lambda \right)=1$.

The equations for each region are

**Region 1:** $0\le \mathit{x}\le 1$

${\mathit{v}}^{\prime}=\frac{\mathit{C}-1}{\mathit{n}}$,

${\mathit{C}}^{\prime}=\frac{\mathrm{vC}-\mathit{x}}{\eta}$.

**Region 2:** $1\le \mathit{x}\le \lambda $

${\mathit{v}}^{\prime}=\frac{\mathit{C}-1}{\mathit{n}}$,

${\mathit{C}}^{\prime}=\frac{\mathrm{vC}-1}{\eta}$.

The interface point $\mathit{x}=1$ is included in both regions. At this point, the solver produces both *left* and *right *solutions, which must be equal to ensure continuity of the solution.

To solve this system of equations in MATLAB, you need to code the equations, boundary conditions, and initial guess before calling the boundary value problem solver `bvp5c`

. You either can include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.

The equations for ${\mathit{v}}^{\prime}\left(\mathit{x}\right)$ and ${\mathit{C}}^{\prime}\left(\mathit{x}\right)$ depend on the region being solved. For multipoint boundary value problems the derivative function must accept a third input argument `region`

, which is used to identify the region where the derivative is being evaluated. The solver numbers the regions from left to right, starting with 1.

Create a function to code the equations with the signature `dydx = f(x,y,region,p)`

, where:

`x`

is the independent variable.`y`

is the dependent variable.`dydx(1)`

gives the equation for ${\mathit{v}}^{\prime}\left(\mathit{x}\right)$, and`dydx(2)`

the equation for ${\mathit{C}}^{\prime}\left(\mathit{x}\right)$.`region`

is the number of the region where the derivative is being computed (in this two-region problem,`region`

is`1`

or`2`

).`p`

is a vector containing the values of the constant parameters $\left[\mathit{n},\text{\hspace{0.17em}}\kappa ,\text{\hspace{0.17em}}\lambda ,\text{\hspace{0.17em}}\eta \right]$.

Use a switch statement to return different equations depending on the region being solved. The function is

function dydx = f(x,y,region,p) % equations being solved n = p(1); eta = p(4); dydx = zeros(2,1); dydx(1) = (y(2) - 1)/n; 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; end end

*Note: All functions are included as local functions at the end of the example.*

Solving two first-order differential equations in two regions requires four boundary conditions. Two of these conditions come from the original problem:

$\mathit{v}\left(0\right)=0$,

$\mathit{C}\left(\lambda \right)-1=0$.

The other two conditions enforce the continuity of the left and right solutions at the interface point $\mathit{x}=1$:

${\mathit{v}}_{\mathit{L}}\left(1\right)-{\mathit{v}}_{\mathit{R}}\left(1\right)=0$,

${\mathit{C}}_{\mathit{L}}\left(1\right)-{\mathit{C}}_{\mathit{R}}\left(1\right)=0$.

For multipoint BVPs, the arguments of the boundary conditions function `YL`

and `YR`

become matrices. In particular, the `k`

th column `YL(:,k)`

is the solution at the left boundary of the `k`

th region. Similarly, `YR(:,k)`

is the solution at the right boundary of the `k`

th region.

In this problem, $\mathit{y}\left(0\right)$ is approximated by `YL(:,1)`

, while $\mathit{y}\left(\lambda \right)$ is approximated by `YR(:,end)`

. Continuity of the solution at $\mathit{x}=1$ requires that `YR(:,1) = YL(:,2)`

.

The function that encodes the residual value of the four boundary conditions is

function res = bc(YL,YR) 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

For multipoint BVPs, the boundary conditions are automatically applied at the beginning and end of the interval of integration. However, you must specify double entries in `xmesh`

for the other interface points. A simple guess that satisfies the boundary conditions is the constant guess `y = [1; 1]`

.

xc = 1; xmesh = [0 0.25 0.5 0.75 xc xc 1.25 1.5 1.75 2]; yinit = [1; 1]; sol = bvpinit(xmesh,yinit);

Define the values of the constant parameters and put them in the vector `p`

. Provide the function to `bvp5c`

with the syntax `@(x,y,r) f(x,y,r,p)`

to provide the vector of parameters.

Calculate the solution for several values of $\kappa $, using each solution as the initial guess for the next. For each value of $\kappa $, calculate the value of the osmolarity $\mathrm{Os}=\frac{1}{\mathit{v}\left(\lambda \right)}$. For each iteration of the loop, compare the computed value with the approximate analytical solution.

lambda = 2; n = 5e-2; for kappa = 2:5 eta = lambda^2/(n*kappa^2); p = [n kappa lambda eta]; sol = bvp5c(@(x,y,r) f(x,y,r,p), @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

2 1.462 1.454 3 1.172 1.164 4 1.078 1.071 5 1.039 1.034

Plot the solution components for $\mathit{v}\left(\mathit{x}\right)$ and $\mathit{C}\left(\mathit{x}\right)$ and a vertical line at the interface point $\mathit{x}=1$. The displayed solution for $\kappa =5$ results from the final iteration of the loop.

plot(sol.x,sol.y(1,:),'--o',sol.x,sol.y(2,:),'--o') line([1 1], [0 2], 'Color', 'k') legend('v(x)','C(x)') title('A Three-Point BVP Solved with bvp5c') xlabel({'x', '\lambda = 2, \kappa = 5'}) ylabel('v(x) and C(x)')

Listed here are the local helper functions that the BVP solver `bvp5c`

calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

function dydx = f(x,y,region,p) % equations being solved n = p(1); eta = p(4); dydx = zeros(2,1); dydx(1) = (y(2) - 1)/n; 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; end end %------------------------------------------- function res = bc(YL,YR) % boundary 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 %-------------------------------------------