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;