% Lorentz

clear all;

sigma = 10;
rho   = 28;
beta  = 8/3;

fun = @(x)( [ sigma*(x(2)-x(1)); rho*x(1) - x(1)*x(3) - x(2);x(1)*x(2)-beta*x(3)]);

p(1,:) = fsolve(fun,[8,8,28])
p(2,:) = fsolve(fun,[-8,-8,28])

p = p';

figure;
hold on;
for i=1:2

  [V,D]=eig(Dapprox(fun,p(:,i)));
  
  EU = V(:,1);
 
  % integrate;
  
  [t,x] = ode45( @(t,x)( -fun(x)/( norm(fun(x)) ) ),[0:0.1: 1000],[p(:,i)+1e-6*EU]);
  
  plot3(x(:,1),x(:,2),x(:,3),'r-');
  
  [t,x] = ode45( @(t,x)( -fun(x)/( norm(fun(x)) ) ),[0:0.1: 1000],[p(:,i)-1e-6*EU]);
  plot3(x(:,1),x(:,2),x(:,3),'r-');
  
end

p = [0;0;0];

[V,D]=eig(Dapprox(fun,p));
  
EU = V(:,1);

% integrate;

[t,x] = ode45( @(t,x)( -fun(x)/( norm(fun(x)) ) ),[0:0.1: 1000],[p+1e-6*EU]);

plot3(x(:,1),x(:,2),x(:,3),'k--');


[t,x] = ode45( @(t,x)( -fun(x)/( norm(fun(x)) ) ),[0:0.1: 1000],[p-1e-6*EU]);

plot3(x(:,1),x(:,2),x(:,3),'k--');