Example 4.8

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

Contents

Settings and Function
ODE Solver
Population Plot
Photon Density Plot

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