% HW3: Henon bifurcation diag for Compu Expt 2.2 p.76
% barnett 10/15/09

N = 300;                                     % # iters (any more and is slow)
x = zeros(2,N+1); b = -0.3;
figure; disp('please wait about 8 seconds...')
for a=0:0.01:2.2
  f = @(x) [a-x(1,:).^2+b*x(2,:);  x(1,:)];  % Henon map w/ params a,b
  x(:,1) = [.5;.5];                          % initial value, better than (0,0)
  for n=1:N, x(:,n+1) = f(x(:,n)); end       % iterate
  plot(a, x(1, end-100:end), 'k.'); hold on; % only plot the final 100 iters
end

% Now, for some fun: want to make it faster? Vectorize by evolving all
% the different cases of a at once! We avoid explicit (thus slow) double loop:
clear
tic;                                               % start the timer
N = 300; b = -0.3;
a = 0:0.01:2.2;
x = [.5;.5] * ones(size(a));
figure;
for n=1:N, x = [a-x(1,:).^2+b*x(2,:); x(1,:)];     % note: a is list
  if n>=N-100, plot(a, x(1,:), 'k.'); hold on; end % add to plot per iteration
end
toc  % elapsed time is now 0.27 secs, a factor 30 faster than other method!
     % you may even decrease a step to 0.001 and it still takes only 0.8 sec
