%% Extend BVP Solution with Extrapolation
% Solve a BVP over an initial interval, and then iteratively extend the interval using each solution as the initial guess for the next interval.
% Consider the equation
% {y^{\prime } }^{\prime } =y.
% As a first-order system, the equation becomes a system of two equations
% {y_1 }^{\prime } =y_2,
% {y_2 }^{\prime } =y_1.
% The equation is initially defined on the interval \left\lbrack 0,3\right\rbrack and is subject to the boundary conditions
% y\left(0\right)=0,
% y\left(3\right)=1.
% The function that encodes the equation as a first-order system is
% function dydx = bvpfun(x,y)
% dydx = [y(2)
%         y(1)];
% end
% Similarly, the function that encodes the boundary conditions is
% function res = bcfun(ya,yb)
% res = [ya(1)
%        yb(1)-1];
% end
% You either can include the required functions as local functions at the end of a file (as done here), or you can save them as separate, named files in a directory on the MATLAB path.
%
%% Initial Guess
% Use an exponential function as the initial guess for the solution. Since the equation has two solution components, write an initial guess function of the form y = guess(x) that returns a vector.
% function y = guess(x)
% y = [exp(x)
%      exp(x)];
% end
% A mesh of five points is sufficient to capture the behavior of the guess function.

xmesh = linspace(0,3,5);
solinit = bvpinit(xmesh,@guess);

%
%% Solve Equation
% Solve the equation in the initial interval \left\lbrack 0,3\right\rbrack and plot the results for y_1.

sol = bvp4c(@bvpfun, @bcfun, solinit);
plot(sol.x(1,:),sol.y(1,:),'-o')

%
%% Extend Interval
% Now, use bvpinit to extend the interval of integration in a loop, solving and plotting each new problem. In each iteration, form the initial guess using the previous solution sol extrapolated to the new interval [0 k]. In each new problem, bvp4c enforces the boundary conditions at the new boundaries [0 k].

hold on
for k = 4:8
    solinit = bvpinit(sol,[0 k]);
    sol = bvp4c(@bvpfun, @bcfun, solinit);
    plot(sol.x(1,:),sol.y(1,:),'-o')
end

% This example shows a simplified version of continuation, a useful technique to solve BVPs by breaking the problem down into smaller intervals or simpler problems. For more examples of this technique, see:
% Solve BVP Using Continuation
% Verify BVP Consistency Using Continuation
%
%% Local Functions
% Listed here are the local helper functions that the BVP solver bvp4c calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

function dydx = bvpfun(x,y) % equation being solved
dydx = [y(2)
        y(1)];
end
%-------------------------------------------
function res = bcfun(ya,yb) % boundary conditions
res = [ya(1)
       yb(1)-1];
end
%-------------------------------------------
function y = guess(x) % guess at solution behavior
y = [exp(x)
     exp(x)];
end
%-------------------------------------------

% Copyright 2018 The MathWorks, Inc.