% visualize & animate phase plane in conservative ODE in potential well.
% barnett 11/9/07

if 0  % .......... pendulum
  P = @(x) -cos(x);          % potential function (eg for pendulum)
  F = @(x) -sin(x);          % force function - make sure F  = -dP/dx
  xl = -5; xh = 10; xs = xl:0.1:xh;            % domain of x
  ys = -4:0.1:4;             % domain of y = x'
  dT = 0.15;
else % .......... cubic half-well (worksheet)
  P = @(x) x - x.^3;          % potential function (eg for pendulum)
  F = @(x) -1 + 3*x.^2;          % force function - make sure F  = -dP/dx
  xl = -3; xh = 3.5; xs = xl:0.05:xh;            % domain of x
  ys = -8:0.2:8;             % domain of y = x'
  dT = 0.03;
end

E = @(x,y) y.^2/2 + P(x);              % total energy function
[xx yy] = meshgrid(xs,ys);
EE = E(xx,yy);
figure; surf(xs,ys,EE); hold on; contour3(xs,ys,EE);
xlabel('x'); ylabel('x'''); zlabel('E(x,x'')'); shading interp; axis vis3d


figure; % animation of flowing particles...
N = 20;
f = @(t,x) [x(2); F(x(1))];     % 1st-order coupled system
x = rand(2,N)*2 - 1; x(1,:) = linspace(xl, xh, N);
for T=0:dT:100
  for n=1:N
    [t, y] = ode45(f, [T T+dT], x(:,n), odeset('abstol', 1e-2));
    x(:,n) = y(end,:)';
  end
  contour(xs,ys, EE); axis tight; hold on; plot(x(1,:), x(2,:), '.');
  xlabel('x'); ylabel('x'''); drawnow; hold off;
end


