%Exercicio Aula 3 
clc; clear all; close all;
%Definição das Variáveis
m=0.1; l=0.2; k=0.015/(0.2^2); g=9.81;
npont=60;
xp = [-3 3;-5 5];  %definição dos limites dos eixos do gráfico e dos estados
x1min = xp(1,1);  %definição de pontos
x1max = xp(1,2);
x2min = xp(2,1);
x2max = xp(2,2);
deltx1=(x1max-x1min)/npont; % Resolução dos estados
deltx2=(x2max-x2min)/npont;

for i=1:2

    x1e(i)=(i-1)*pi;
    a=[0 1;-g*cos(x1e(i))/l -k/m];  %montando a matriz para calcular os autovalores
    disp('auto-valor da matriz a para i igual a:')
    i  %exibir os pontos de equilíbrio
    eig(a) %calcular os autovalores de A
end
%resolução das EDO´S

t=[0:0.001:10]; %tempo de resoluçao da EDO
y0=[0.5 1];    %condição inicial para Y
[T,Y]=ode23(@rigid,t,y0);

%plotar V e a solução

figure (1)
q = [-2:deltx1:2];
w = [-5:deltx2:5];
[q,w] = meshgrid(q,w);
z=(1-cos(q)).*g/l+w.^2/2;         
contour3 (q,w,z,100,'b');
hold on
plot(Y(1,1),Y(1,2),'r*',Y(:,1),Y(:,2),'r')
xlabel('x1');
hold off
ylabel('x2');
zlabel('Função de Energia');

%plotar V2 no tempo

p12=0.5*k/m;
p11=k/m*p12;
p22=1;
p21=p12;
P=[p11 p12;p21 p22];
for i=1:size(Y,1)
V(i)=0.5*Y(i,:)*P*Y(i,:)'+(g/l)*(1-cos(Y(i,1)));
end
figure(2)
plot(T,V)
xlabel('Tempo');
ylabel('Energia');

%determinar a positividade de V2

syms p11 p12 p22  %comando para declarar variável
P=[p11 p12; p12 p22];  %montamos matriz de variáveis
eig(P)       %calcular autovalores