Listings for nonlinear dynamics purpose

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)])

%{
File: test_phasePortrait.m
Author: Jonathan C. Teixeira
Email: jonta.teixeira@gmail.com
Description: Generate phase portrait of ode system
%}

% as usual
clc; clear; close all;

ode = 'nonlinear';    % ode selections: 'linear' or 'nonlinear'

switch lower(ode)
case 'linear'
    % homogeneous linear system
    % name of variables
    varname = {'x', 'y'};
    portrait_lim = {[-5,-5],[5,5]};
    num_curves = 100;
    [x0, y0] = deal(0,0);

    a   = [-1, +1, -1, +1];
    b   = [-5, -5, +5, +5];

    for i=1:numel(a)
      % define ode system
      dy = @(x) [a(i)*x(1); b(i)*(x(1) + x(2))];

      % plot phase portrait
      phasePortrait(dy,portrait_lim{1},portrait_lim{2},num_curves)
      hold on
      plot(x0,y0,'sr');
      xlabel(varname{1}); ylabel(varname{2}); grid on
      print('-dpng', '-r300', ['phase_portrait_',ode,'_',num2str(i)]);
    end

case 'nonlinear'
    %% nonlinear system
    % parameter definition
    tau = 60;
    k0  = 1.7e5;
    G   = 90.3;
    E_R = 6307;
    Ca0 = 1.5;
    T0  = linspace(250,350,12);

    % name of variables
    varname = {'Concentration', 'Temperature'};
    portrait_lim = {[0,200],[3,550]};
    num_curves = 300;

    for i=1:numel(T0)
        % define ode system
        dy = @(x) [(Ca0 - x(1))/tau - k0*exp(-E_R/x(2))*x(1); ...
                   (T0(i) - x(2))/tau + G*k0*exp(-E_R/x(2))*x(1)];

        % critical points
        opts = optimoptions('fsolve','Display','off','MaxFunctionEvaluations',3000, 'functiontolerance',1e-6);
        x = [fsolve(dy,[0;0], opts), ...
                fsolve(dy,[0;500], opts), ...
                fsolve(dy,[1;380], opts)]
        j = [sqrt(sum(dy(x(:,1)).^2))<1e-3; ...
             sqrt(sum(dy(x(:,2)).^2))<1e-3; ...
             sqrt(sum(dy(x(:,3)).^2))<1e-3]
        disp(i);
        y0 = x(2,j); x0 = x(1,j);


        % plot phase portrait
        phasePortrait(dy,portrait_lim{1},portrait_lim{2},num_curves, false)
        hold on
        plot(x0,y0,'sr','markerfacecolor','r');
        xlabel(varname{1}); ylabel(varname{2}); grid on;
        title(['T0 = ',num2str(T0(i))])
        print('-dpng', '-r300', ['phase_portrait_',ode,'_',num2str(i)]);
    end
otherwise
  error('** select valid model')
end