ch2ex9

function ch2ex9
d = 0.1;
yd = 1 - (1/2)*d^2 - (1/6)*d^4;
erryd = (19/180)*d^6;
fprintf('At d = %g, the error in y(d) is about %5.1e.\n',d,erryd);
[y,x] = ode45(@ode,[yd 0],d);
fprintf('Total collapse occurs at x = %g.\n',x(end));
% Augment the arrays with y = 1 at x = 0 for the plot
% and plot with the original independent variable x.
y = [1; y];
x = [0; x];
%plot(x,y)
plot(x,y,'-k')
%print -depsc ch2fig9
%========================================================
function dxdy = ode(y,x)
dxdy = - sqrt(3*y^3/(2*(1 - y^3)));