%%
%  Author: Peter Polcz (ppolcz@gmail.com) 
%  Created on 2023. October 12. (2023a)
%

syms t x1 x2 real
x = [x1;x2];

A = sym('a%d%d',[2,2]);
b = sym('b%d',[2,1]);

f = x .* (A*x + b);

J = jacobian(f,x);

% Normalized vector field
Norm_f = norm(f);

f1 = f(1) / Norm_f;
f2 = f(2) / Norm_f;

f_fh = matlabFunction(f,'Vars',{x,A,b});
f1_fh = matlabFunction(f1,'Vars',{x1,x2,A,b});
f2_fh = matlabFunction(f2,'Vars',{x1,x2,A,b});

J_fh = matlabFunction(J,'Vars',{x,A,b});

clear x A b J f1 f2

%%

params = linspace(-0.1,0.1,51);

for p = params
    
    A = [
           0.000    p
           1.000    0.000
        ];
    
    b = [
           1.600
          -1.200
        ];
    
    
    Eq0 = [0;0];
    Eq1 = -A\b;
    
    LaTeX = sprintf([ '$'...
        '\\left\\{\\matrix{ ' ...
        '    \\dot x = x (%g %+g y) ' ...
        '    \\cr' ...
        '    \\dot y = y (%g %+g x) ' ...
        '}\\right. ~~~ (p = %g)' ...
        '$'],b(1),A(1,2),b(2),A(2,1),p);
    
    EqX = [-A(1,1)\b(1) ; 0];
    EqY = [0 ; -A(2,2)\b(2)];
    
    Eq = [ Eq0 Eq1 EqX EqY ];
    
    % Remove equilibrium points, which do not belong the the positive quadrant
    Eq(:,any(Eq < 0,1)) = [];
    
    Lim = 1.2*max(Eq(:));
    % Lim = 8;
    lspace = linspace(0,Lim,51);
    [x1,x2] = meshgrid(lspace);
    
    f1 = f1_fh(x1,x2,A,b);
    f2 = f2_fh(x1,x2,A,b);
    
    [t,x] = ode45(@(t,x) f_fh(x,A,b),[0 100],[0.1;0.1]);

    fig = figure(312);
    delete(fig.Children);
    Tl = tiledlayout(1,2,'TileSpacing','compact','Padding','compact');
    ax = nexttile;
    hold on, grid on, box on;
    
    quiver(x1,x2,f1,f2,'DisplayName','Vector field')
    
    for i = 1:size(Eq,2)
        plot(Eq(1,i),Eq(2,i),'.','MarkerSize',25)
    end
  
    plot(x(1,1),x(1,2),'r.');
    plot(x(:,1),x(:,2),'r');

    axis equal
    title(LaTeX,'Interpreter','latex','FontSize',14)
    ax.XLim = [0,Lim];
    ax.YLim = [0,Lim];

    ax = nexttile;
    hold on, grid on, box on;
    
    Tr = linspace(-3,3,101);
    Pl = plot(Tr,Tr.^2/4);
    for i = 1:size(Eq,2)
        J = J_fh(Eq(:,i),A,b);
        d = det(J);
        t = trace(J);
    
        plot(t,d,'.','MarkerSize',25)
    end
    XLim = ax.XLim;
    xgrid = linspace(XLim(1),XLim(2),100000);
    YLim = ax.YLim;
    plot(XLim,[0 0],'Color',Pl.Color);
    plot([0 0],YLim,'Color',Pl.Color);
    plot(xgrid,xgrid.^2/4,'Color',Pl.Color);
    ax.XLim = XLim;
    ax.YLim = YLim;

    drawnow
    keyboard
end