function Alinearized_example

figure(8); clf
xfixed = [1;2]
for n=1:40
%    x0 = (rand(2,1)-.5)*.1 + xfixed
    x0 = (rand(2,1)-.5)*10 + xfixed;
    if n>30, x0 = ginput(1); end
    % NUMERICAL INTEGRATION EXAMPLE:
    %S = odeset('AbsTol',1e-6,'RelTol',1e-6);
    S = odeset('AbsTol',1e-3,'RelTol',1e-3);
    [tout,xout] = ode45(@Alin_ode,[0 200],x0,S);
    %[tout,xout] = ode45(@Alin_ode,[0 1000],x0,S);
    plot(xout(:,1),xout(:,2),'b'); hold on
    plot(xout(end,1),xout(end,2),'r.')
    axis([-10 10 -10 10]*.5 + [xfixed(1)+[0 0], xfixed(2)+[0 0]])
    
%    plot(xout(end-100:end,1),xout(end-100:end,2),'r-')
end

function dX = Alin_ode(t,X,plot_nc)
a = 0; b = 1; c = -4; d = -1;  % A: damped oscillation
%a = 0; b = 1; c = -4; d = -.05; % B: light damping
%a = 0; b = 1; c = -4; d = -4; % C: critically damped
%a = 0; b = 1; c = -4; d = -6; % D: over-damped
a = 0; b = 1; c = -4; d = 0; % E: no damping?
%a = -4; b = 0; c = 0; d = -2;  % F: stable, uncoupled
%a = -2; b = 0; c = 0; d = -2; % G: STAR node

x = X(1);
y = X(2);
f = a*(x-1)+(x-1)^2 + b*sin(y-2);
g = (1/9)*(x-1)^3 + c*(x-1) + (d/2)*(sin((y-2))+(y-2));
dX = [f;g];

% Below, linearized system:
%dX = [a, b; c, d]*(X - [1;2]);
