brach.m

function brach
%BRACH    Brachistochrone illustration.
%   Computes and plots approximate brachistochrone by optimization,
%   using fminsearch, and exact brachistochrone, using fzero.

bx = 1; g = 9.81;
byvals = linspace(0.2,2,10);
Nvals = [4 8];
for i = 1:2
   N = Nvals(i);
   subplot(1,2,i)
   for k = 1:length(byvals)

      % Approximate brachistochrone.
      by = byvals(k);
      dy = by/N; dx = bx/N;
      yinit = [dy:dy:by-dy];
      y = fminsearch(@Btime,yinit);

      plot([0:dx:bx],-[0 y by],'ro-')
      hold on

      % True brachistochrone.
      tzero = @(theta)(by*theta - by*sin(theta) + bx*cos(theta) - bx);
      tstar = fzero(tzero,pi);
      R = by/(1-cos(tstar));
      thetavals = linspace(0,tstar,100);
      xcoord = R*(thetavals-sin(thetavals));
      ycoord = R*(1-cos(thetavals));
      plot(xcoord,-ycoord,'b--','Linewidth',2)

   end
   title(sprintf('N = %1.0f',N),'FontWeight','normal','FontSize',12)
   xlim([0,bx]), axis off
end
hold off

  function T = Btime(y)
  %BTIME    Travel time for a particle.
  %   Piecewise linear path with equispaced y between (0,0) and (bx,by).

  yvals = [0 y by];        % End points do not vary.
  N = length(y)+1; d = bx/N;
  T = sum(2*sqrt( d^2 + (diff(yvals)).^2 )./( sqrt(2*g*yvals(2:end)) + ...
          sqrt(2*g*yvals(1:end-1) )));
  end

end