clear all;

% Constant Mass
opt = odeset ("RelTol", 1e-4, "AbsTol", 1e-5, ...
              "Jacobian", [98, 198; -99, -199], "Mass", eye(2), ...
              "MStateDependence", "none");
[t, y] = ode15s (@(t, y)[98, 198; -99, -199] * (y - [1; 0]), ...
                 [0, 5], [2; 0], opt);
y1xct = @(t) 2 * exp (-t) - exp (-100 * t) + 1;
y2xct = @(t) - exp (-t) + exp (-100 * t);
assert ([y1xct(t), y2xct(t)], y, 1e-3);
figure; hold on;
plot (t, y(:, 1), '.'); plot (t, y(:, 2), '.');

clear all;

% Jacobian as fun
opt = odeset ("RelTol", 1e-4, "AbsTol", 1e-5, ...
              "Jacobian", [98, 198; -99, -199]);
[t, y] = ode15s (@(t, y)[98, 198; -99, -199] * (y - [1; 0]), ...
                 [0, 5], [2; 0], opt);
y1xct = @(t) 2 * exp (-t) - exp (-100 * t) + 1;
y2xct = @(t) - exp (-t) + exp (-100 * t);
assert ([y1xct(t), y2xct(t)], y, 1e-3);
figure; hold on;
plot (t, y(:, 1), '.'); plot (t, y(:, 2), '.');

clear all;

% Jacobian as const sparse
opt = odeset ("RelTol", 1e-4, "AbsTol", 1e-5, ...
              "Jacobian", sparse([98, 198; -99, -199]));
[t, y] = ode15s (@(t, y)[98, 198; -99, -199] * (y - [1; 0]), ...
                 [0, 5], [2; 0], opt);
y1xct = @(t) 2 * exp (-t) - exp (-100 * t) + 1;
y2xct = @(t) - exp (-t) + exp (-100 * t);
assert ([y1xct(t), y2xct(t)], y, 1e-3);
figure; hold on;
plot (t, y(:, 1), '.'); plot (t, y(:, 2), '.');

% Jacobian as const => ode15s.m BIST
test ode15s

