ex37cch2.m

these codes are for solving OED with matlab

function ex37cch2
% Emden problem
global n
options = odeset('Events',@events,'RelTol',1e-8,'AbsTol',1e-10);
n = 3;
xend   = 10;
xinit  = 0.1
yinit  = 1 - (1/6)*xinit^2 + (n/120)*xinit^4;
ypinit = -(1/3)*xinit + (4*n/120)*xinit^3;
erry   = ((5*n-8*n*n)/(3*5040))*xinit^6
erryp  = (6*(5*n-8*n*n))/(3*5040)*xinit^5
[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)')
title(['Emden problem for xinit = 'num2str(xinit)])
xinit2  = .05
yinit2  = 1 - (1/6)*xinit2^2 + (n/120)*xinit2^4;
ypinit2 = -(1/3)*xinit2 + (4*n/120)*xinit2^3;
erry2   = ((5*n-8*n*n)/(3*5040))*xinit2^6
erryp2  = (6*(5*n-8*n*n))/(3*5040)*xinit2^5
[x2,y2,xe2,ye2,ie2] = ode15s(@f,[xinit2 xend],[yinit2; ypinit2],options);
if isempty(ie2)
    fprintf(' y(x) was never 0.')
else
    fprintf(' y(x)=0 for x = %g.\n',x2(end))
end
figure
plot(x2,y2(:,1),x2,y2(:,2))
legend('y_1(t)','y_2(t)')
title(['Emden problem for xinit = 'num2str(xinit2)])
%===================================================================
function dydt = f(x,y)
global n
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = -(2/x)*y(2) - (y(1))^n;

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