function phasePortrait(fun, lb, ub, n, randomize)
%Create a phase portrait from a two-dimensional system of ode
%
% SYNOPSIS:
%   phasePortrait(fun, lb, ub, n, randomize)
%
% PARAMETERS:
%   fun         - Function handle or function name with the system of ode 
%   lb          - lower bound of phase portrait
%   ub          - upper bound of phase portrait
%   n           - number of intial condition of system
%   randomize   - logical variable to randomize the initial points to integrate 
%                 the system. OPTIONAL.  Default value == true
%
% RETURNS:
%   phase portrait plot
%
% EXAMPLE:
%

if nargin<4 
  n = 100;
  randomize = true;
elseif nargin < 5
  randomize = true;
end

num = 20;
dx = abs(ub-lb)/num;

[x1,x2]=meshgrid(lb(1):dx(1):ub(1),lb(2):dx(2):ub(2));

v = zeros(numel(x1(:)),2);
for i=1:numel(x1(:))
   v(i,:) = fun([x1(i), x2(i)]);
end

figure()
disp('vectorial field')
hold on
quiver(x1(:),x2(:),v(:,1),v(:,2));
if randomize
  u=lb(1)+(ub(1)-lb(1))*rand(n,1);
  v=lb(2)+(ub(2)-lb(2))*rand(n,1);
else 
  n = floor(n/4+1);
  u = linspace(lb(1),ub(1),n); 
  u = [repmat(u,1,2),repmat(lb(1),1,n), repmat(ub(1),1,n)]; 
  v = linspace(lb(2),ub(2),n);
  v = [repmat(lb(2),1,n), repmat(ub(2),1,n), repmat(v,1,2)]; 
end

funcao = @(t,x) fun([x(1),x(2)]);
disp('tracing paths ...')
for i=1:numel(u)
    [~,sol]=ode45(funcao,[0,1000],[u(i),v(i)]);
    plot(sol(:,1),sol(:,2),'k');
end
hold off

xlim([lb(1) ub(1)])
ylim([lb(2) ub(2)])
