%% Example 4.8
% A Stiff System (Byrne and Hindmarsh)
% Here we use ode23s

%% Settings and Function
close all;
clear a b g r s u f J t y options;

a = 1.5e-18;
b = 2.5e-6;
g = 2.1e-6;
r = .6;
s = .18;
u = .016;

% function
f = @(t,y) [-y(1)*(a*y(2)+b)+g;y(2)*(r*y(1)-s)+u*(1+y(1))];

% Jacobian
J = @(t,y) [-(a*y(2)+b) -a*y(1); r*y(2)+u r*y(1)-s];

%% ODE Solver
% Here we turn off Matlab's warning.  This problem leads to a poorly
% conditioned linear system.
warning off all
options = odeset('RelTol',1e-4,'AbsTol',1e-4*ones(2,1),'Jacobian',J);
[t,y] = ode23s(f,[0 7e5],[-1 0]',options);
warning on all

%% Population Plot
figure(1)
plot(t,y(:,1))
xlabel('Time')
ylabel('Population Inversion')
grid on

%% Photon Density Plot
figure(2)
semilogy(t,y(:,2))
xlabel('Time')
ylabel('Photon Density')
grid on