ex19ch3.m

these codes are for solving OED with matlab

function sol = ex19ch3
global alpha epsilon
alpha = 0.5;
epsilon = 1;
solinit = bvpinit(linspace(0,1,20),@guess,pi);
sol = bvp4c(@odes,@bcs,solinit);
plot(sol.x,sol.y(1,:),[0 1],[0 0])
title(['Whirling string with \alpha = ',num2str(alpha),...
       ', \epsilon = ',num2str(epsilon)])
fprintf('Computed omega = %g.\n',sol.parameters);
S = deval(sol,0.5);
fprintf('Computed y(0.5) is %g.\n',S(1));
S = deval(sol,1);
fprintf('Computed y(1.0) is %g.\n',S(1));
%==========================================
function v = guess(x)
v = zeros(4,1);
v(1) = cos(pi*x);
v(2) = -pi*sin(pi*x);
v(3) = x;
v(4) = 1;

function dydx = odes(x,y,omega)
global alpha epsilon
dydx = zeros(4,1);
temp = sqrt(1 + epsilon^2 *y(1)^2);
dydx(1) = y(2);
dydx(2) = - omega^2*y(1)*(alpha^2 + (1 - alpha^2)/(y(4)*temp));
dydx(3) = 1/temp;
dydx(4) = 0;

function res = bcs(ya,yb,omega)
global alpha epsilon
res = [ya(1) - 1
       ya(2)
       ya(3)
       yb(2)
       yb(4) - (1 -(1 - alpha^2)*yb(3))/alpha^2];