clear all 
close all 
clc 

% initial conditions span

% Simulacao nao-linear
x10 = -10:2:10;
x20 = -10:2:10;

% Plano de fase 1: pto equilibrio (0,0)
% x10 = -4-pi:2:4-pi;
% x20 = -10:2:10;

% Plano de fase 2: pto equilibrio (pi,0)
% x10 = -pi:1:pi;
% x20 = -10:0.5:10;

for i=1:length(x10)
    for j=1:length(x20)
        x0 = [x10(i);x20(j)];
        
        [t,x] = ode45(@plano_P1,[0:0.05:20],x0); % time simulation
        figure(1)
        plot(t,x(:,1),t,x(:,2));
        %axis([0 20 -30 30])
        hold on
        
        figure(2)
        plot(x(:,1),x(:,2),'-b','LineWidth',1)
        %axis([-20 20 -20 20])
        hold on 
        
    end
end
figure(2)  
% Plano de fase nao-linear

% Plano de fase 1
axis([-10 10 -10 10])

%Plano de fase 2
axis([-pi pi -10 10])

function xdot = plano_nao_linear(t,x)
g = 9.8;
ell = 1;
m = 1;
b = 1;

x1 = x(1);
x2 = x(2);

x1dot = x2;
x2dot = -(g/ell)*sin(x1)-(b/m)*x2;
xdot = [x1dot;x2dot];
end

% Ponto de equilibrio (0,0)
function xdot = plano_P1(t,x)
g = 9.8;
ell = 1;
m = 1;
b = 1;

x1 = x(1);
x2 = x(2);

x1dot = x2;
x2dot = -(g/ell)*x1-(b/m)*x2;

xdot = [x1dot;x2dot];
end

% Ponto de equilibrio (pi,0)
function xdot = plano_P2(t,x)
g = 9.8;
ell = 1;
m = 1;
b = 1;

x1 = x(1);
x2 = x(2);

x1dot = x2;
x2dot = (g/ell)*x1-(b/m)*x2;

xdot = [x1dot;x2dot];
end

