// Equation definition deff('[r]=chemres(t,y,yd)',[ 'r(1)=-0.04*y(1)+1d4*y(2)*y(3)-yd(1);'; 'r(2)=0.04*y(1)-1d4*y(2)*y(3)-3d7*y(2)*y(2)-yd(2);' 'r(3)=y(1)+y(2)+y(3)-1;']) deff('[p]=chemad(t,y,p)','p=p+diag([1 1 0])') deff('[p]=chemjac(t,y,yd)',[ 'p=[-0.04, 1.d4*y(3) , 1.d4*y(2);'; ' 0.04, -1d4*y(3)-6d7*y(2) , -1d4*y(2);'; ' 1 , 1 , 1 ]']) // Integration y0=[1;0;0]; yd0=[-0.04;0.04;0]; t=[1.d-5:0.02:.4 0.41:.1:4 40 400 4000 40000 4d5 4d6 4d7 4d8 4d9 4d10]; rtol = 1d-4;atol=[1.d-6;1.d-10;1.d-6]; y=impl(y0,yd0,0,t,rtol,atol,chemres,chemad,chemjac);
// visualization 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)