function void = stableBVP(void)
% First example of Chapter 1.  
set(gcf,'DefaultAxesColorOrder',[0 0 0; 0 0 0; 0 0 0])
% Now compute the equilibrium solution as a BVP.

  for infinity = 95:100,
  plot([0,infinity],[0,0])
  axis([0 infinity -3 8])
  hold on
  solinit = bvpinit(linspace(0,infinity,10),[3 0]);
  sol = bvp4c(@pendulum,@pendulumbc,solinit);

  plot(sol.x,sol.y(1,:),':ko','MarkerSize',2,'MarkerFaceColor','k')

  disp(['With infinity = ', num2str(infinity)])
  fprintf('Initial slope of equilibrium solution is computed to be %e.\n',sol.y(2,1));
  disp('Hit return to continue')
  pause
end

%=============================================
function yp = pendulum(t,y)
yp = [y(2); -sin(y(1))];

function res = pendulumbc(ya,yb)
res = [    ya(1)
	    yb(1) - pi ]; 


