function [t,f,fa] = sampleproblem3(tmax,threshold);
%[t,f,fa] = sampleproblem3(tmax,threshold)
%
%SOLVE the linear coupled ode pair:
%[df1dt;df2dt] = [1 -2;0 -5][f1;f2]
%with initial condition: [f1_0;f2_0] = [10;10]
%which has the analytical solution
%f1(t) = 15exp(-t)-5exp(-5t)
%f2(t) = 10exp(-5t)
%
%stop integrating when f2 reaches value "threshold"
%
%Show mastry of plotting options
%
%VARIABLES:
%tmax = integration time
%t = column vector of integration timepoints
%f = numerical solution at timepoints:   f(:,1)  = f1(t);  f(:,2) = f2(t)
%fa = analytical solution at timepoints: fa(:,1) = f1a(t); f(:,2) = f2a(t)
%threshold = f2 value that stops integration early
%
%EXAMPLE
%[t,f,fa] = sampleproblem3(5,1);


param.threshold = threshold;
options = odeset('Events',@eventfunction);
f0 = [10;10];
[t,f] = ode45(@derivfunction,[0;tmax],f0,options,param);
fa = [15*exp(-t)-5*exp(-5*t), 10*exp(-5*t)]; %fa not plotted in this example

%NEW MATERIAL: "help plot" gives more options on format strings
figure(1);
plot(t,f(:,1),'r^',t,f(:,2),'--b');
legend('f1','f2');
ymax = max(max(f));
ymin = 0;
xmax = max(t);
xmin = 0;
xlim([xmin,xmax]);
ylim([ymin,ymax]);
help plot
return;

function derivs = derivfunction(t,f,param)
derivs = [-1 2;0 -5]*f;
return;

function [value,isterminal,direction] = eventfunction(t,f,param);
    value = f(2) - param.threshold;
    isterminal = 1;
    direction = 0;
return;