%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];
x1min = xp(1,1);
x1max = xp(1,2);
x2min = xp(2,1);
x2max = xp(2,2);
deltx1=(x1max-x1min)/npont;
deltx2=(x2max-x2min)/npont;

for i=1:2

    x1e(i)=(i-1)*pi;
    a=[0 1;-g*cos(x1e(i))/l -k/m];
    disp('auto-valor da matriz a para i igual a:')
    i
    eig(a)
end
%resolução das EDO´S

t=[0:0.001:10];
y0=[0.5 1];
[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('funcao 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)

%determinar a positividade de V2

syms p11 p12 p22
P=[p11 p12; p12 p22];
eig(P)



    