%
% 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