%
%   BESS_EX1.M  File to Plot Solution to a Specific BVP
%
%   This goal of this file is to solve the following BVP in terms of Bessel functions.
%      y'' + 9xy = 0   with   y(1) = -10   and   y(2) = 0 
%
%   The general solution is given in terms of J and Y Bessel functions of order 1/3
%   as follows (see class notes):
%      y(x) = sqrt(x)*(C1*J(1/3)(2x^1.5) + C2*x*Y(1/3)(2x^1.5))
%   and the C1 & C2 coefficients for the given BCs are given by the 
%   solution to the following system of equations (xo = 1 & xf = 2):
%      sqrt(xo)*J(1/3)(2*xo^1.5)*C1  +  sqrt(xo)*Y(1/3)(2*xo^1.5)*C2  =  y(xo)
%      sqrt(xf)*J(1/3)(2*xf^1.5)*C1  +  sqrt(xf)*Y(1/3)(2*xf^1.5)*C2  =  y(xf))
%
%   Just for fun, a numerical solution using the Shooting Method was added for 
%   comparison purposes (calls bvp2sh.m).
%
%   File prepared by J. R. White, UMass-Lowell (Nov. 2003)
%

%
%   getting started
      clear all,  close all
%
%   ANALYTICAL SOLUTION using Bessel Functions
%
%   evaluate coefficients 
      xo = 1;   xf = 2;   yo = -10;   yf = 0;
      A = [sqrt(xo)*besselj(1/3,2*xo^1.5)   sqrt(xo)*bessely(1/3,2*xo^1.5);
           sqrt(xf)*besselj(1/3,2*xf^1.5)   sqrt(xf)*bessely(1/3,2*xf^1.5)];
      b = [yo yf]';
      C = A\b;
%
%   now evaluate the unique solution over the domain of interest
      xa = linspace(xo,xf,51);
      ya = sqrt(xa).*(C(1)*besselj(1/3,2*xa.^1.5) + C(2)*bessely(1/3,2*xa.^1.5));
%
%   NUMERICAL SOLUTION using the Shooting Method
%
%   solution via Shooting Method (numerical approx)
      zbc = [1 0 yo;            % left BC  --> y(xo) = yo
             1 0 yf];           % right BC --> y(xf) = yf
      tol = 1e-3;  options = odeset('RelTol',tol);   % set tolerance for ODE soln
      [xs,zs] = bvp2sh('bess_ex1a',[xo xf],zbc,options);
      ys = zs(:,1);    % store solution for plotting
%
%   plot both solutions
      figure(1)
      plot(xa,ya,'r-',xs,ys,'bo','LineWidth',2),grid
      title('Bess\_Ex1:  Solution to specific BVP using Bessel Functions')
      xlabel('x values'),ylabel('y(x)')
      legend('Analytical Soln','Numerical Soln')
%
%   end of problem