clear all
close all
clc

% global mi

% Parametros do gráfico

    n = 10000;     % Número de amostragem
    escala = .7;    % Tamaho das setas
    npont=50;       % Número de pontos em cada eixo
    
    
    % Área onde será plotado os vetores
    xp = [-10 10;-10 10];
    x1min = xp(1,1);
    x1max = xp(1,2);
    x2min = xp(2,1);
    x2max = xp(2,2);
    
    deltx1=(x1max-x1min)/npont;
    deltx2=(x2max-x2min)/npont;
    
    compr=sqrt(deltx1^2+deltx2^2)/2;

    i=0;
    j=0;
    
% Desenhar o Campo de Vetores

for x1=x1min:deltx1:x1max  % Varia de x1min a x1max com passo deltx1
    i=i+1;

for x2=x2min:deltx2:x2max % Varia de x2min a x2max com passo deltx2
    x = [x1 x2]; t=0;
    
% Calcula o xponto para t = 0 e x = [x1 x2]
xponto = equacao(t,x); 
m = xponto(2)/(xponto(1)+1e-8);     %corresponde ao sentido da f(x)
dx1 =compr/sqrt(1+m^2);
j=j+1;
dx2 = m*dx1;
if x2 < 0
dx2 = -dx2;
dx1 = -dx1;
end
dx(j,i) = dx1*.5;
dy(j,i) = dx2*.5;
hold on;
end
j=0;
end
% Mapa vetorial
[x,y] = meshgrid(x1min:deltx1:x1max,x2min:deltx2:x2max);

% Plota os vetores de coordenadas (dx/2) e (dy/2) nos pontos x e y
quiver(x,y,dx*.5,dy*.5,escala);  
axis([x1min x1max x2min+deltx2 x2max+deltx2]);

% Define o vetor tempo (t)
tinicial = 0;
tfinal = 100;
t=linspace(tinicial,tfinal,n); % cria o vetor tempo (t) com 100000 pontos, de 0 à 100 segundos.

% CÁLCULO DAS CONDIÇÕES INICIAIS QUE SERÃO UTILIZADAS NA ODE45.

N=5;   % N e M definem a quantidade de condições iniciais que serão consideradas.
M=5;   % O total de soluções (em vermelho) que serão plotadas é dado pelo 
        % produto N*M.
   
R = 1; % utilizado para calcular as N*M condições iniciais

for I = 1:N+1,
for J = 1:M+1,
x0(R,1) = x1min+((x1max-x1min)/N)*(I-1);  % cria as condições iniciais referentes à variável de estado x1
x0(R,2) = x2min+((x2max-x2min)/M)*(J-1);  % cria as condições iniciais referentes à variável de estado x2
R=R+1;
end
end

% Obtem as trajetórias para a matriz de condiçoes iniciais x0 
hold on;
for i = 1:N*M,
[tempo,x] = ode45('equacao',t,x0(i,:));
plot(x(:,1),x(:,2),'r'); % plota todas as soluções de acordo com cada condição inicial                
end