% Boundary Value Problems
% To demonstrate the use of finite difference methods for solving one-dimensional boundary value problems, we will solve the following problem:
% \frac{d^2 y}{dx^2 }+y=1
% with the two boundary conditions:
% \begin{array}{l}
% y\left(0\right)=1\\
% y\left(\frac{\pi }{2}\right)=0
% \end{array}
% The solution to this differential equation is
% y\left(x\right)=1-\sin \left(x\right)
% So how do we solve this using finite differences? Here we divide the region 0<x<\pi /2 into several intervals. For our purposes we will choose five intervals, giving us a spacing between mesh points of h=\pi /10. We will refer to these points as x=x_i, where the first point is x_0  and the last point is  x_5. This allows us to write x_{i+1} =x_i +h. Similarly we can refer to the y-values at each of the mesh points as y_i =y\left(x_i \right). So our goal is to find the y-values at each of the mesh points. In this example, there are six of them.
% Two of the y-values can be found from the boundary conditions. That is, y_0 =1 and y_4 =0. So now we just need to find the four internal values. To do this we approximate the derivative in terms of the internal y-values. At some point halfway between x_i and x_{i+1} we can approximate the derivative of y as:
% {\left|\frac{\textrm{dy}}{\textrm{dx}}\right|}_{i+\frac{1}{2}} \approx \frac{y_{i+1} -y_i }{h}
% Similarly, the derivative at some point half-way between xi-1and xi can be approximated as:
% {\left|\frac{\textrm{dy}}{\textrm{dx}}\right|}_{i-\frac{1}{2}} \approx \frac{y_i -y_{i-1} }{h}
% Then we can approximate the second derivative at point xi as:
% {\left|\frac{d^2 \mathrm{y}}{\mathrm{d}{\mathrm{x}}^2 }\right|}_i \approx \frac{{\left|\frac{\textrm{dy}}{\textrm{dx}}\right|}_{i+\frac{1}{2}} -{\left|\frac{\textrm{dy}}{\textrm{dx}}\right|}_{i-\frac{1}{2}} }{h}
% or
% {\left|\frac{d^2 \mathrm{y}}{\mathrm{d}{\mathrm{x}}^2 }\right|}_i \approx \frac{y_{i+1} -2y_i +y_{i-1} }{h^2 }
% Now we can just plug this into our model equation, giving us:
% \frac{y_{i+1} -2y_i +y_{i-1} }{h^2 }+y_i =1
% This difference equation can then be written at the four internal mesh points, giving us four algebraic equations for the four unknowns. The methods used to complete this process depend on the tool to be used.
% Note that if you have to solve a different differential equation, you have to re-derive this formula.
% To solve our problem using MatLab, we construct a matrix equation and then solve it using the built-in solver. To facilitate this, we write the difference equation as:
%  y_{i+1} -\left(2-h^2 \right)y_i +y_{i-1} =h^2
% By writing this at the three internal mesh points, we obtain:
%  \begin{array}{l}
% y_0 -\left(2-h^2 \right)y_1 +y_2 =h^2 \\
% y_1 -\left(2-h^2 \right)y_2 +y_3 =h^2 \\
% y_2 -\left(2-h^2 \right)y_3 +y_4 =h^2
% \end{array}
% Since we know the boundary values, we can write:
% \begin{array}{l}
% y_0 =1\\
% y_0 -\left(2-h^2 \right)y_1 +y_2 =h^2 \\
% y_1 -\left(2-h^2 \right)y_2 +y_3 =h^2 \\
% y_2 -\left(2-h^2 \right)y_3 +y_4 =h^2 \\
% y_4 =0
% \end{array}
% Now we would like to write this as a matrix equation:
%  \begin{array}{l}
% \textrm{Ax}=B\\
% A=\left\lbrack \begin{array}{ccccc}
% 1 & 0 & 0 & 0 & 0\\
% 1 & -\left(2-h^2 \right) & 1 & 0 & 0\\
% 0 & 1 & -\left(2-h^2 \right) & 1 & 0\\
% 0 & 0 & 1 & -\left(2-h^2 \right) & 1\\
% 0 & 0 & 0 & 0 & 1
% \end{array}\right\rbrack \\
% B=\left\lbrack \begin{array}{c}
% 1\\
% h^2 \\
% h^2 \\
% h^2 \\
% 0
% \end{array}\right\rbrack
% \end{array}
% This equation is then easily solved using Matlab. Note that if we want to increase the number of divisions, we just change h, increase the size of the A matrix and adjust A and B appropriately.
% My solution is as follows:

npoints=5;
b=pi/2;
h=b/npoints;
A=[1 0 0 0 0; ...
    1 -(2-h^2) 1 0 0;...
    0 1 -(2-h^2) 1 0;...
    0 0 1 -(2-h^2) 1;...
    0 0 0 0 1]
B=[1; h^2; h^2; h^2; 0]
y=A\B
x=linspace(0,b,npoints);
plot(x,y,'LineWidth',2)
xlabel('x')
ylabel('y')
set(gca,'FontSize',20)
exact=1-sin(x);
hold on
plot(x,exact,'o')
hold off

% To generalize this to more points (and therefore obtain more accuracy), we just have to find clever ways to produce the proper A and B. Here is one example

npoints=50;
b=pi/2;
h=b/npoints;
vn=ones(npoints,1);
vnm1=ones(npoints-1,1);
A=-(2-h^2)*diag(vn)+diag(vnm1,1)+diag(vnm1,-1);
A(1,1)=1;A(1,2)=0;
A(npoints,npoints-1)=0;A(npoints,npoints)=1;
B=ones(npoints,1)*h^2;
B(1)=1; B(npoints)=0;
y=A\B;
x=linspace(0,b,npoints);
plot(x,y,'LineWidth',2)
xlabel('x')
ylabel('y')
set(gca,'FontSize',20)
exact=1-sin(x);
hold on
plot(x,exact,'o')
hold off

% A Built-In Matlab Routine(bvp4c)
% Matlab includes the routine bvp4c to help solve boundary value problems. This routine uses finite differences to solve a problem, adjusting the mesh in order to obtain desired accuracy. The user provides a routine (odefun) that defines a system of first order equations and a second function (bcfun) that defines a set of boundary conditions at each end of the interval of interest. The user also must provide locations of the initial guesses for the dependent variables, allowing for the use of nonlinear equations. Matlab then returns solutions for all the dependent variable. The solution, using this routine, for the model problem in this documents is given here:

xlow=0;
xhigh=pi/2;
%use 100 mesh points
x=linspace(xlow,xhigh,100);
%use guess of y=1 and dy/dx=-1 at all points
solinit = bvpinit(x,[1 -1]);
sol = bvp4c(@bvp4ode,@bvp4bc,solinit);
%evaluate solution
Sxint = deval(sol,x);
plot(x,Sxint(1,:),'LineWidth',2)
xlabel('x')
ylabel('solution y')
set(gca,'FontSize',20)
exact=1-sin(x);
hold on
plot(x,exact,'o')
hold off
% -----------------------------------------------

%

function dydx = bvp4ode(x,y,lambda)
%this function defines system of two first order
%equations representing our original second
%order equation
%dy/dx=z
%dz/dx=1-y
%so y(1) represents y from the original equation
%and y(2) represents dy/dx
dydx = [  y(2)     1-y(1) ];
end
% -----------------------------------------------
function res = bvp4bc(ya,yb,lambda)
%This defines y(0)=1 and y(pi/2)=0
res = [  ya(1)-1   yb(1) ];
end

%  In this script, the bvp4c function defines the boundary conditions, the bvp4ode function defines the first order ordinary differential equations that result from the differential equation, and the bvp4 function organizes the calculation. Within bvp4, the solinit variable defines the mesh using the linspace function, the sol variable carries the solution after running the bvp4c function, and Sxint evaluates the result at various values of x.
% © James Blanchard, University of Wisconsin-Madison, blanchard@engr.wisc.edu