%Exercicio Aula 3 
clc; clear all; close all;
%Definição das Variáveis
c=2*10^-3; r=5*10^3;l=5*10^3; e=sqrt(l/c);
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;
  % x1e(1)=(0.0063);
  %  a=[-(1/c)*(17.16-2*103.79*(x1e)+229*(x1e.^2)-905.2*(x1e.^3)+418.6*(x1e.^4)) 1/c;-1/l -r/l];  %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=linspace(0,5,30); %tempo de resoluçao da EDO
y0=[0.5 -0.3];    %condição inicial para Y
[T,Y]=ode45(@rigid_diodo_tunel,t,y0);

%/plotar V e a solução

figure (1)
q = [-.6:deltx1:.6];
w = [-.6:deltx2:.6];
[q,w] = meshgrid(q,w);

eps=0.01;
hy1=17.76*q-103.79*q.^2+229.62*q.^3-226.31*q.^4+83.72*q.^5;

z=(1/2)*c*( q.^2 + (eps.*hy1 + w).^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');

figure(2)

plot(0,Y(1,1),'r*',T,Y(:,1));

figure(3)

plot(Y(1,1),Y(1,2),'r*',Y(:,1),Y(:,2),'r')
grid