function f = equilibrium_points(type,x0)

%type [string] is the type of equilibrium

T = 3;
dt = 0.05;

switch num2str(type)
    case {'stable focus'}
        disp('stable focus')
        a = -1;
        b = 2;
        c = -2;
        d = -1;
    case {'unstable focus'}
        disp('unstable focus')
        a = 1;
        b = 2;
        c = -2;
        d = 1;
    case {'center'}
        disp('center')
        a = 1;
        b = 2;
        c = -2;
        d = -1;
    case {'stable node'}
        disp('stable node')
        a = -1;
        b = -0.5;
        c = -0.5;
        d = -1;
    case {'unstable node'}
        disp('unstable node')
        a = 1;
        b = -0.5;
        c = -0.5;
        d = 1;
    case {'saddle point'}
        disp('saddle point')
        a = 1;
        b = -1.5;
        c = -1.5;
        d = 1;
    case {'stable line'}
        disp('stable line')
        a = -1;
        b = -1.0;
        c = -1.0;
        d = -1;
    case {'unstable line'}
        disp('unstable line')
        a = 1;
        b = -1.0;
        c = -1.0;
        d = 1;
    case {'stable star'}
        disp('stable star')
        a = -1;
        b = 0;
        c = 0;
        d = -1;
    case {'unstable star'}
        disp('unstable star')
        a = 1;
        b = 0;
        c = 0;
        d = 1;
    case {'degenerate node'}
        disp('degenerate node')
        a = 1;
        b = 1;
        c = 0;
        d = 1;
    case {'trivial'}
        disp('trivial')
        a = 0;
        b = 1;
        c = 0;
        d = 0;
    otherwise
        error('unknown choice of equilibrium')
end
A = [a b; c d]
[V D] = eig(A);
d = diag(D)

B = [0; 0];
C = [0 0];
D = 0;

t = 0:dt:T;
Nt = length(t);
u = zeros(1,Nt);
SYS = ss(A,B,C,D);
[y,t,x] = lsim(SYS,u,t,x0);


if isreal(d(1))
    plot([0 V(1,1)],[0 V(2,1)],'r','LineWidth',2)
    hold on
    plot([0 V(1,2)],[0 V(2,2)],'r','LineWidth',2)
else
    hold on
end

if strcmp(type(1:2),'un')
    xM = 10;
    yM = 10;
else
    xM = 4;
    yM = 4;
end
set(gca,'XLim',[-xM xM],'Ylim',[-yM yM], ...
    'DataAspectRatio',[1 1 1]);

for i = 1:Nt
    plot(x(1:i,1),x(1:i,2),'-')
    pause(0.01)
end
title(num2str(type))