%%
%  Author: Peter Polcz (ppolcz@gmail.com) 
%  Created on 2023. December 07. (2023a)
%
% [1] Kocaoglu, Karal and Guzelis (2013). Analysis of chaotic dynamics of Chua's circuit
% with lncosh nonlinearity. 2013 8th International Conference on Electrical and
% Electronics Engineering (ELECO), IEEE, 2013.
% 
% Chua circuit bifurcation analysis. Periodic solutions, chaotic solutions: single scroll,
% double scroll.

% In this parameter setup, the physical parameters R, C1, C2, L are given indirectly.
alpha = 9;
beta = 100/7;
gamma = 0.016;
m0 = -8/7;
m1 = -5/7;

% `abs(x+1) - abs(x-1)` is approximated by the smooth `lncosh(x+1) - lncosh(x-1)` 
lncosh = @(lambda,x) lambda\log(cosh(lambda*x));

lambda = 2.88; % 1.5, 1.8, 1.9, 1.95, 2.88, 3.3740955 4
h = @(x) m1*x + 2\(m0-m1)*(lncosh(lambda,x+1) - lncosh(lambda,x-1));

f = @(t,x) [
    alpha * ( x(2) - x(1) - h(x(1)) )
    x(1) - x(2) + x(3)
    -beta*x(2) - gamma*x(3)
    ];

% First solve the ODE, then, consider the last state vector as the initial condition for
% the next solution. (In this way, the periodic solutions can be simply detected.)
[t,x] = ode89(f,[0,100],[1;1;1]/10);
[t,x] = ode89(f,[0,1000],x(end,:)');

fig = figure(2);
Tl = tiledlayout(3,2,'TileSpacing','tight','Padding','compact');

% Plot the solution in the phase diagram.
ax = nexttile([3,1]);
plot3(x(:,1),x(:,2),x(:,3))

xlabel x
ylabel y
zlabel z
grid on, box on
axis vis3d

% Plot the x1, x2, and x3 functions
ax1 = nexttile; grid on, box on, hold on
plot(t,x(:,1))

ax2 = nexttile; grid on, box on, hold on
plot(t,x(:,2))

ax3 = nexttile; grid on, box on, hold on
plot(t,x(:,3))

Link = linkprop([ax1,ax2,ax3],'XLim');
xlim([0,200])

