% ### LINode45EX.m ###       01.26.16 CB
% Numerically integrate a general 2nd order linear autonomous system (w/
% const. coefficients)
% x' = a*x + b*y
% y' = c*x + d*y

clear
% -----------------------------------------------------
% User input (Note: All paramters are stored in a structure)
P.y0(1) = 1.0;   % initial value for x 
P.y0(2) = 1;   % initial value for y
P.A= [-3.9 3;      % matrix A to contain coefficients A= [a b
      -2 1];     %                                      c d]

% Integration limits
P.t0 = 0.0;   % Start value
P.tf = 10.0;   % Finish value
P.dt = 0.01;  % time step
% ----------------------------------------------------------------------
% +++
% determine some basic derived quantities 
p= P.A(1,1)+ P.A(2,2);  % Tr(A)
q= P.A(1,1)* P.A(2,2)-P.A(1,2)* P.A(2,1); % det(A)
disp(['Tr(A)= ' num2str(p),' and det(A)= ',num2str(q)]); 
eigV1= [0.5*(p+sqrt(p^2-4*q)) 0.5*(p-sqrt(p^2-4*q))];     % calc. eigenvalues directly
eigV2= eig(P.A);    % calculate via Matlab's built-in routine
disp(['eigenvalues= ' num2str(eigV1(1)),' and ',num2str(eigV1(2))]);

% +++
% use built-in ode45 to solve
[t y] = ode45('LINfunction', [P.t0:P.dt:P.tf],P.y0,[],P);

% ------------------------------------------------------
% visualize
% NOTE (re variable naming): x=y(1) and y=y(2)
figure(1); clf;
plot(t,y(:,1)); hold on; grid on; 
xlabel('t');    ylabel('x(t)')
% Phase plane 
figure(2); clf;
plot(y(:,1), y(:,2)); hold on; grid on; 
xlabel('x(t)');    ylabel('y(t)')
% "solution space" 
figure(3); clf;
plot(p,q,'rx','MarkerSize',9,'LineWidth',3); hold on; grid on;
if (abs(p)<1), pSpan= linspace(-1,1,100);
else    pSpan= linspace(-1.5*p-0.1,1.5*p+0.1,100); end
qSpan= pSpan.^2/4;
plot(pSpan,qSpan,'k-','LineWidth',2); %ylim([-max(qSpan) max(qSpan)])
plot(pSpan,zeros(numel(pSpan),1),'b--','LineWidth',2);
xlabel('Tr(A)');    ylabel('det(A)')