%% Solve BVP with Unknown Parameter
% This example shows how to use bvp4c to solve a boundary value problem with an unknown parameter.
% Mathieu's equation is defined on the interval \left\lbrack 0,\pi \right\rbrack by
% {y^{\prime } }^{\prime } +\left(\lambda -2q\;\cos \left(2x\right)\right)y=0.
% When the parameter q=5, the boundary conditions are
% y^{\prime } \left(0\right)=0,
% y^{\prime } \left(\pi \right)=0.
% However, this only determines y\left(x\right) up to a constant multiple, so a third condition is required to specify a particular solution,
% y\left(0\right)=1.
% 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 bvp4c. You can either 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 Equation
% Create a function to code the equations. This function should have the signature dydx = mat4ode(x,y,lambda), where:
% x is the independent variable.
% y is the dependent variable.
% lambda is an unknown parameter representing an eigenvalue.
% You can write Mathieu's equation as a first-order system using the substitutions y_1 =y and y_2 =y^{\prime },
% {y_1 }^{\prime } =y_2,
% {y_2 }^{\prime } =-\left(\lambda -2q\;\cos \left(2x\right)\right)y_1.
% The corresponding function is then
% function dydx = mat4ode(x,y,lambda) % equation being solved
% dydx = [y(2)
%       -(lambda - 2*q*cos(2*x))*y(1)];
% end
% Note: All functions are included as local functions at the end of the example.
%
%% Code Boundary Conditions
% Now, write a function that returns the residual value of the boundary conditions at the boundary points. This function should have the signature res = mat4bc(ya,yb,lambda), where:
% ya is the value of the boundary condition at the beginning of the interval \left\lbrack a,b\right\rbrack.
% yb is the value of the boundary condition at the end of the interval \left\lbrack a,b\right\rbrack.
% lambda is an unknown parameter representing an eigenvalue.
% This problem has three boundary conditions in the interval \left\lbrack 0,\pi \right\rbrack. To calculate the residual values, you need to put the boundary conditions into the form g\left(x,y\right)=0. In this form the boundary conditions are
% y^{\prime } \left(0\right)=0,
% y^{\prime } \left(\pi \right)=0,
% y\left(0\right)-1=0.
% The corresponding function is then
% function res = mat4bc(ya,yb,lambda) % boundary conditions
% res = [ya(2)
%        yb(2)
%        ya(1)-1];
% end
%
%% Create Initial Guess
% Lastly, create an initial guess of the solution. You must provide an initial guess for both solution components y_1 =y\left(x\right) and y_2 =y^{\prime } \left(x\right), as well as the unknown parameter \lambda.
% For this problem, a cosine function makes for a good initial guess since it satisfies the three boundary conditions. Code the initial guess for y using a function that returns the guess for y_1 and y_2.
% function yinit = mat4init(x) % initial guess function
% yinit = [cos(4*x)
%         -4*sin(4*x)];
% end
% Call bvpinit using a mesh of 10 points in the interval \left\lbrack 0,\pi \right\rbrack, the initial guess function, and a guess of 15 for the value of \lambda.

lambda = 15;
solinit = bvpinit(linspace(0,pi,10),@mat4init,lambda);

%
%% Solve Equation
% Call bvp4c with the ODE function, boundary condition function, and initial guess.

sol = bvp4c(@mat4ode, @mat4bc, solinit);

%
%% Value of Parameter
% Print the value of the unknown parameter \lambda found by bvp4c. This value is the fourth eigenvalue (q=5) of Mathieu's equation.

fprintf('Fourth eigenvalue is approximately %7.3f.\n',...
            sol.parameters)

%
%% Plot Solution
% Use deval to evaluate the solution computed by bvp4c at 100 points in the interval \left\lbrack 0,\pi \right\rbrack.

xint = linspace(0,pi);
Sxint = deval(sol,xint);

% Plot both solution components. The plot shows the eigenfunction (and its derivative) associated with the fourth eigenvalue \lambda_4 =17\ldotp 097.

plot(xint,Sxint)
axis([0 pi -4 4])
title('Eigenfunction of Mathieu''s Equation.') 
xlabel('x')
ylabel('y')
legend('y','y''')

%
%% 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 = mat4ode(x,y,lambda) % equation being solved
q = 5;
dydx = [y(2)
      -(lambda - 2*q*cos(2*x))*y(1)];
end
%-------------------------------------------
function res = mat4bc(ya,yb,lambda) % boundary conditions
res = [ya(2)
       yb(2)
       ya(1)-1];
end
%-------------------------------------------
function yinit = mat4init(x) % initial guess function
yinit = [cos(4*x)
        -4*sin(4*x)];
end
%-------------------------------------------

% Copyright 2018 The MathWorks, Inc.