function TimeSeries = chua(L,R0,C2,G,Ga,Gb,C1,E,x0,y0,z0,dataset_size,step_size)
% Syntax: TimeSeries=chua(L,R0,C2,G,Ga,Gb,C1,E,x0,y0,z0,dataset_size,step_size)
% ______________________________________
%
% Chua's oscillator is described by a  set  of  three
% ordinary differential equations called Chua's equations:
%
%      dI3      R0      1
%      --- =  - -- I3 - - V2
%      dt       L       L
%
%      dV2    1       G
%      --- = -- I3 - -- (V2 - V1)
%      dt    C2      C2
%
%      dV1    G              1
%      --- = -- (V2 - V1) - -- f(V1)
%      dt    C1             C1
%
% where
%
%   f(V1) = Gb V1 + - (Ga - Gb)(|V1 + E| - |V1 - E|)
%
% A solution of these equations (I3,V2,V1)(t) starting from an
% initial state (I3,V2,V1)(0) is called a trajectory of Chua's
% oscillator.
%
% This function uses a RungeKutta integration method optimised for the chua
% paradigm
%
% Reference:
% ABC - Adventures in Bifurication & Chaos ... Prof M.P Kennedy 1993
%
%
% James McEvoy, Tom Murray
% 
% University e-mail: 99375940@student.ucc.ie
% Lifetime e-mail: sacevoy@eircom.net
% Homepage: http://www.sacevoy.com
%
% 2 Nov 2002



% Models Initial Variables
%-------------------------

TimeSeries = [x0, y0, z0]'; % models initial conditions

% Optimized Runge-Kutta Variables
%--------------------------------

h = step_size*G/C2;
h2 = (h)*(.5);
h6 = (h)/(6);

anor = Ga/G;
bnor = Gb/G;
bnorplus1 = bnor + 1;
alpha = C2/C1;
beta = C2/(L*G*G);
gammaloc = (R0*C2)/(L*G);

bh = beta*h;
bh2 = beta*h2;
ch = gammaloc*h;
ch2 = gammaloc*h2;
omch2 = 1 - ch2;

k1 = [0 0 0]';
k2 = [0 0 0]';
k3 = [0 0 0]';
k4 = [0 0 0]';
M = [0 0 0]';

% Calculate Time Series
%----------------------

M(1) = TimeSeries(3)/E;
M(2) = TimeSeries(2)/E;
M(3) = TimeSeries(1)/(E*G);

for i=1:dataset_size
    % Runge Kutta
    % Round One
    k1(1) = alpha*(M(2) - bnorplus1*M(1) - (.5)*(anor - bnor)*(abs(M(1) + 1) - abs(M(1) - 1)));
    k1(2) = M(1) - M(2) + M(3);
    k1(3) = -beta*M(2) - gammaloc*M(3);
    % Round Two
    temp = M(1) + h2*k1(1);
    k2(1) = alpha*(M(2) + h2*k1(2) - bnorplus1*temp - (.5)*(anor - bnor)*(abs(temp + 1) - abs(temp - 1)));
    k2(2) = k1(2) + h2*(k1(1) - k1(2) + k1(3));
    k2(3) = omch2*k1(3) - bh2*k1(2);
    % Round Three
    temp = M(1) + h2*k2(1);
    k3(1) = alpha*(M(2) + h2*k2(2) - bnorplus1*temp - (.5)*(anor - bnor)*(abs(temp + 1) - abs(temp - 1)));
    k3(2) = k1(2) + h2*(k2(1) - k2(2) + k2(3));
    k3(3) = k1(3) - bh2*k2(2) - ch2*k2(3);
    % Round Four
    temp = M(1) + h*k3(1);
    k4(1) = alpha*(M(2) + h*k3(2) - bnorplus1*temp - (.5)*(anor - bnor)*(abs(temp + 1) - abs(temp - 1)));
    k4(2) = k1(2) + h*(k3(1) - k3(2) + k3(3));
    k4(3) = k1(3) - bh*k3(2) - ch*k3(3);
    
    M = M + (k1 + 2*k2 + 2*k3 + k4)*(h6);
    
    TimeSeries(3,i+1) = E*M(1);
    TimeSeries(2,i+1) = E*M(2); 
    TimeSeries(1,i+1) = (E*G)*M(3);
    
    i=i+1;
end