# Solve BVP with Multiple Boundary Conditions

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.

### Code Equations

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.*

### Code Boundary Conditions

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

### Form Initial Guess

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);

### Solve Equation

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 Solution

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)')

### Local Functions

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 %-------------------------------------------