Chemical process (Stiff)

-->// Equations definition
-->deff('[yd]=chem(t,y)',[
-->    'yd(1)=-0.04*y(1) + 1d4*y(2)*y(3);';
-->    'yd(3)= 3d7*y(2)*y(2);';
-->    'yd(2)= -yd(1) - yd(3);'])
-->// Integration
-->t=[1.d-5:0.02:.4 0.41:.1:4 40 400 4000 40000 4d5 4d6 4d7 4d8 4d9 4d10];
-->rtol=1.d-4;atol=[1.d-6;1.d-10;1.d-6];
-->y=ode([1;0;0],0,t,rtol,atol,chem);
-->// Visualisation

-->rect=[1.d-5,-0.1,1d11,1.1];
-->plot2d1("oln",t',(diag([1 10000 1])*y)',(1:3),"111",' y1@10000 y2@y3',rect)
-->// Add surface condition
-->nt=prod(size(t));
-->deff('[y]=surf(t,x)','y=[x(1)-1.e-4;x(3)-1.e-2]')
-->// First root
-->[y,rd,w,iw]=ode('root',[1;0;0],0,t,rtol,atol,chem,2,surf);rd
 rd  =
!   0.2640459    2. !
-->while rd<>[] then
-->  [nw,ny]=size(y);
-->  k=find(rd(1)>t(1:nt-1)&rd(1)<t(2:nt));
-->  // Visualisation
-->  write(%io(2),[rd(1);y(:,ny)]','(''t='',e10.3,'' y='',3(e10.3,'',''))')
-->  plot2d1("oln",rd(1)',(diag([1 10000 1])*y(:,ny))',[-3,-3,-3],"000");
-->  // Next root
-->  [y,rd,w,iw]=ode('root',[1;0;0],rd(1),t(k+1:nt),rtol,atol,chem,2,surf,w,iw);
-->end
t= 0.264E+00 y= 0.990E+00, 0.347E-04, 0.100E-01,
t= 0.207E+08 y= 0.100E-03, 0.400E-09, 0.100E+01,