%% Compare bvp4c and bvp5c Solvers
% Solve a BVP at a crude error tolerance with two different solvers and compare the results.
% Consider the second-order ODE
% {y^{\prime } }^{\prime } +\frac{2}{x}y^{\prime } +\frac{1}{x^4 }y=0.
% The equation is defined on the interval \left\lbrack \frac{1}{3\pi },1\right\rbrack subject to the boundary conditions
% y\left(\frac{1}{3\pi }\right)=0,
% y\left(1\right)=\sin \left(1\right).
% To solve this equation in MATLAB®, you need to write a function that represents the equation as a system of first-order equations, write a function for the boundary conditions, set some option values, and create an initial guess. Then the BVP solver uses these four inputs to solve the equation.
%
%% Code Equation
% With the substitutions y_1 =y and y_2 =y^{\prime }, you can rewrite the ODE as a system of first-order equations
% {y_1 }^{\prime } =y_2,
% {y_2 }^{\prime } =-\frac{2}{x}y_2 -\frac{1}{x^4 }y_1.
% The corresponding function is
% function dydx = bvpfcn(x,y)
% dydx = [y(2)
%        -2*y(2)/x - y(1)/x^4];
% end
% Note: All functions are included at the end of the example as local functions.
%
%% Code Boundary Conditions
% The boundary condition function requires that the boundary conditions are in the form g\left(y\left(a\right),y\left(b\right)\right)=0. In this form, the boundary conditions are
% y\left(\frac{1}{3\pi }\right)=0,
% y\left(1\right)-\sin \left(1\right)=0.
% The corresponding function is
% function res = bcfcn(ya,yb)
% res = [ya(1)
%        yb(1)-sin(1)];
% end
%
%% Set Options
% Use bvpset to turn on the display of solver statistics, and specify crude error tolerances to highlight the difference in error control between the solvers. Also, for efficiency, specify the analytical Jacobian
% J=\frac{\partial f_i }{\partial y}=\left\lbrack \begin{array}{cc}
% \frac{\partial f_1 }{\partial y_1 } & \frac{\partial f_1 }{\partial y_2 }\\
% \frac{\partial f_2 }{\partial y_1 } & \frac{\partial f_2 }{\partial y_2 }
% \end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
% \;\;\;\;0 & \;\;\;\;1\\
% -\frac{1}{x^4 } & -\frac{2}{x}
% \end{array}\right\rbrack.
% The corresponding function that returns the value of the Jacobian is
% function dfdy = jac(x,y)
% dfdy = [0      1
%        -1/x^4 -2/x];
% end

opts = bvpset('FJacobian',@jac,'RelTol',0.1,'AbsTol',0.1,'Stats','on');

%
%% Create Initial Guess
% Use bvpinit to create an initial guess of the solution. Specify a constant function as the initial guess with an initial mesh of 10 points in the interval \left\lbrack 1/3\pi ,1\right\rbrack.

xmesh = linspace(1/(3*pi), 1, 10);
solinit = bvpinit(xmesh, [1; 1]);

%
%% Solve Equation
% Solve the equation with both bvp4c and bvp5c.

sol4c = bvp4c(@bvpfcn, @bcfcn, solinit, opts);

%

sol5c = bvp5c(@bvpfcn, @bcfcn, solinit, opts);

%
%% Plot Results
% Plot the results of the two calculations for y_1 with the analytic solution for comparison. The analytic solution is
% y_1 =\sin \left(\frac{1}{x}\right),
% y_2 =-\frac{1}{x^2 }\cos \left(\frac{1}{x}\right).

xplot = linspace(1/(3*pi),1,200);
yplot = [sin(1./xplot); -cos(1./xplot)./xplot.^2];

plot(xplot,yplot(1,:),'k',sol4c.x,sol4c.y(1,:),'r*',sol5c.x,sol5c.y(1,:),'bo')
title('Comparison of BVP Solvers with Crude Error Tolerance')
legend('True','BVP4C','BVP5C')
xlabel('x')
ylabel('solution y')

% The plot confirms that bvp5c directly controls the true error in the calculation, while bvp4c controls it only indirectly. At more stringent error tolerances, this difference between the solvers is not as apparent.
%
%% Local Functions
% Listed here are the local functions that the BVP solvers use to solve the problem.

function dydx = bvpfcn(x,y) % equation to solve
dydx = [y(2)
       -2*y(2)/x - y(1)/x^4];
end
%---------------------------------
function res = bcfcn(ya,yb) % boundary conditions
res = [ya(1)
       yb(1)-sin(1)];
end
%---------------------------------
function dfdy = jac(x,y) % analytical jacobian for f
dfdy = [0      1
       -1/x^4 -2/x];
end
%---------------------------------

% Copyright 2018 The MathWorks, Inc.