ex37bch2.m

these codes are for solving OED with matlab

function ex37bch2
% Emden problem with luck
global n
options = odeset('Events',@events,'RelTol',1e-8,'AbsTol',1e-10);
n = 3;
xend   = 10;
xinit  = 0;
yinit  = 1;
ypinit = 0;
[x,y,xe,ye,ie] = ode15s(@f,[xinit xend],[yinit; ypinit],options);
if isempty(ie)
    fprintf(' y(x) was never 0.')
else
    fprintf(' y(x)=0 for x = %g.\n',x(end))
end
plot(x,y(:,1),x,y(:,2))
legend('y_1(t)','y_2(t)')
%===================================================================
function dydt = f(x,y)
global n
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = -1/3;
if x > 0
    dydt(2) = -(2/x)*y(2) - (y(1))^n;
end

function [value,isterminal,direction] = events(x,y)
value = y(1);
isterminal = 1;
direction  = 0;