%% testing ode15s with a stiff system
clear *
%% works well if Jacobian passed as function handle
opt = odeset ("RelTol", 1e-8, "AbsTol", 1e-6,'Jacobian',@(t,y)[98, 198;-99, -199],'Stats','on');
%% does not return results of Jacobian passed as matrix
%opt = odeset ("RelTol", 1e-8, "AbsTol", 1e-6,'Jacobian',[98, 198;-99, -199],'Stats','on');
f = @(t,y)[98, 198;-99, -199]*(y-[1;0]);
y_exact = @(t)2*exp(-t)-exp(-100*t)+1;

disp('results with Jacobian')
[t15sJ, y15sJ] = ode15s (f ,[0, 5],[2;0], opt);  y15sJ = y15sJ(:,1);

opt = odeset ("RelTol", 1e-8, "AbsTol", 1e-6,'Stats','on');
disp('results without Jacobian')
[t15s, y15s] = ode15s (f ,[0, 5],[2;0], opt);  y15s = y15s(:,1);

figure(1)
plot(t15s,y_exact(t15s),t15sJ,y15sJ,t15s,y15s)
xlabel('time t'); ylabel('solution y(t)');
legend('exact','with J','no J')
title('solution')

figure(2)
plot(t15sJ,y15sJ-y_exact(t15sJ),t15s,y15s-y_exact(t15s))
xlabel('time t');
legend('with J','no J')
title('difference to exact solution')
