%
%   BESS_EX2.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.
%      xy'' + 2y' - y = 0   with   y(3) = 2   and   y'(6) = -2 
%
%   The general solution is given in terms of I1 and K1 Bessel functions 
%   as follows (see class notes):
%      y(x) = (C1*I1(2*sqrt(x)) + C2*K1(2*sqrt(x)))/sqrt(x)
%   and the C1 & C2 coefficients for the given BCs are given by the 
%   solution to the system of equations derived in class (this was rather 
%   complicated  --  see class notes):
%
%   As a check, 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 = 3;   xf = 6;   yo = 2;   ypf = -2;
      a11 = besseli(1,2*sqrt(xo))/sqrt(xo);
      a12 = besselk(1,2*sqrt(xo))/sqrt(xo);
      a21 =  besseli(0,2*sqrt(xf))/xf - besseli(1,2*sqrt(xf))/xf^1.5;
      a22 = -besselk(0,2*sqrt(xf))/xf - besselk(1,2*sqrt(xf))/xf^1.5;
      A = [ a11   a12;  a21   a22];
      b = [yo ypf]';
      C = A\b;
%
%   now evaluate the unique solution over the domain of interest
      xa = linspace(xo,xf,31);
      ya = (C(1)*besseli(1,2*sqrt(xa)) + C(2)*besselk(1,2*sqrt(xa)))./sqrt(xa);
%
%   NUMERICAL SOLUTION using the Shooting Method
%
%   solution via Shooting Method (numerical approx)
      zbc = [1 0 yo;            % left BC  --> y(xo) = yo
             0 1 ypf];          % right BC --> y'(xf) = ypf
      tol = 1e-3;  options = odeset('RelTol',tol);   % set tolerance for ODE soln
      [xs,zs] = bvp2sh('bess_ex2a',[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\_Ex2:  Solution to specific BVP using Bessel Functions')
      xlabel('x values'),ylabel('y(x)')
      legend('Analytical Soln','Numerical Soln')
%
%   end of problem