 
//Copyright INRIA
 
 
//                 DASRT
 
//
 
//-----------------------------------------------------------------------
 
// First problem.
 
// The initial value problem is..
 
//   DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y,  Y(1) = 1,  10.LE. T0.LE. 6
 
// The solution is  Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3
 
// The two root functions are..
 
//   G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT)  (with root at T = 2.5),
 
//   G2 = LOG(Y) - 2.2491  (with roots at T = 2.47 and 2.53)
 
//-----------------------------------------------------------------------
 
Eps=1.e-3
 Eps  =
 
    0.001  
 
y0=1;t=2:6;t0=1;y0d=3;
 
info=list([],0,[],[],[],0,0);
 
atol=1.d-6;rtol=0;ng=2;
 
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
 
if abs(nn(1)-2.47) > Eps then bugmes();quit;end
 
y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
 
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
 
if abs(nn(1)-2.5)  > Eps then bugmes();quit;end
 
y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
 
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
 
if abs(nn(1)-2.53) > Eps then bugmes();quit;end
 
 
deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2*log(y)+8)/t-5)*y')
 
deff('[rts]=gr1(t,y,yd)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
 
 
y0=1;t=2:6;t0=1;y0d=3;
 
info=list([],0,[],[],[],0,0);
 
atol=1.d-6;rtol=0;ng=2;
 
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
 
if abs(nn(1)-2.47)  > Eps then bugmes();quit;end
 
y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
 
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
 
if abs(nn(1)-2.5)  > Eps then bugmes();quit;end
 
y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
 
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
 
if abs(nn(1)-2.53)  > Eps then bugmes();quit;end
 
 
//Copyright INRIA
 
 
//-----------------------------------------------------------------------
 
// Second problem (Van Der Pol oscillator).
 
// The initial value problem is..
 
//   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
 
//   Y1(0) = 2,  Y2(0) = 0,  00.LE. T0.LE. 200
 
//   Y1PRIME(0) = 0, Y2PRIME(0) = -2
 
// The root function is  G = Y1.
 
// An analytic solution is not known, but the zeros of Y1 are known
 
// to 15 figures for purposes of checking the accuracy.
 
//-----------------------------------------------------------------------
 
Eps=1.e-5
 Eps  =
 
    0.00001  
 
rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
 
t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
 
info=list([],0,[],[],[],0,0);
 
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
 
if abs(nn(1)-81.163512)  > Eps then bugmes();quit;end
 
 
deff('[delta,ires]=res2(t,y,ydot)',...
'ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]')
 
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,'jac2',ng,'gr2',info);
 
deff('J=jac2(t,y,ydot,c)','y1=y(1);y2=y(2);J=[c,-1;200*y1*y2+1,c-100*(1-y1*y1)]')
 
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info);
 
deff('s=gr2(t,y,yd)','s=y(1)')
 
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info);
 
 
//           Hot Restart
 
 
[yy,nn,hot]=dasrt([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
 
t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
 
[yy,nn,hot]=dasrt([y01,y0d1],t01,t,atol,rtol,'res2','jac2',ng,'gr2',info,hot);
 
if abs(nn(1)-162.57763) > Eps then bugmes();quit;end
 
//Copyright INRIA
 
 
files=G_make(['/tmp/res22.o'],'res22.dll');
 
link(files,'res22');
shared archive loaded
Link done
 
files=G_make(['/tmp/jac22.o'],'jac22.dll');
 
link(files,'jac22');
shared archive loaded
Link done
 
files=G_make(['/tmp/gr22.o'],'gr22.dll');
 
link(files,'gr22');
shared archive loaded
Link done
 
//
 
rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
 
//
 
t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
 
//
 
info=list([],0,[],[],[],0,0);
 
//Calling the routines by dasrt
 
//
 
[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
 
// hot restart (uncomment)
 
//[yy,nn,hot]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
 
//t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
 
//[yy,nn,hot]=dasrt([y01,y0d1],t01,t,atol,rtol,'res22','jac22',ng,'gr22',info,hot);
 
 
//Copyright INRIA
 
//DASSL
 
// PROBLEM 1:   LINEAR DIFFERENTIAL/ALGEBRAIC SYSTEM
 
//
 
//X1DOT + 10.0*X1 = 0
 
//X1 + X2 = 1
 
//X1(0) = 1.0, X2(0) = 0.0
 
//
 
t=1:10;t0=0;y0=[1;0];y0d=[-10;0];
 
info=list([],0,[],[],[],0,0);
 
// with default atol and rtol parameters
 
//    Calling Scilab functions
 
deff('[r,ires]=dres1(t,y,ydot)','r=[ydot(1,:)+10*y(1,:);y(2,:)+y(1,:)-1];ires=0')
 
deff('pd=djac1(t,y,ydot,cj)','pd=[cj+10,0;1,1]')
 
//   scilab function, without jacobian
 
yy0=dassl([y0,y0d],t0,t,dres1,info);
 
if norm(dres1(yy0(1,:),yy0(2:3,:),yy0(4:5,:)),1)>1.d-6 then bugmes();quit;end
 
//   scilab functions, with jacobian
 
yy1=dassl([y0,y0d],t0,t,dres1,djac1,info);
 
if norm(dres1(yy1(1,:),yy1(2:3,:),yy1(4:5,:)),1)>1.d-6 then bugmes();quit;end
 
// fortran routine dres1 in dir. routines/default, without jocabian
 
yy2=dassl([y0,y0d],t0,t,'dres1',info);   //=yy0
 
if norm(dres1(yy2(1,:),yy2(2:3,:),yy2(4:5,:)),1)>1.4d-6 then bugmes();quit;end
 
// fortran routines dres1 and djac1 in dir. routines/default, with jacobian
 
yy3=dassl([y0,y0d],t0,t,'dres1','djac1',info);  //=yy1
 
if norm(dres1(yy3(1,:),yy3(2:3,:),yy3(4:5,:)),1)>1.d-6 then bugmes();quit;end
 
//if norm(yy3-yy1,1)  > Eps then bugmes();quit;end
 
yy3bis=dassl([y0,y0d],t0,t,'dres1',djac1,info);
 
if norm(dres1(yy3bis(1,:),yy3bis(2:3,:),yy3bis(4:5,:)),1)>1.d-6 then bugmes();quit;end
 
// call fortran dres1 and scilab's djac1
 
yy3ter=dassl([y0,y0d],t0,t,dres1,'djac1',info);
 
if norm(dres1(yy3ter(1,:),yy3ter(2:3,:),yy3ter(4:5,:)),1)>1.d-6 then bugmes();quit;end
 
//
 
// with specific atol and rtol parameters
 
atol=1.d-6;rtol=0;
 
//   scilab function, without jacobian
 
yy4=dassl([y0,y0d],t0,t,atol,rtol,dres1,info);
 
if norm(dres1(yy4(1,:),yy4(2:3,:),yy4(4:5,:)),1)>1.d-5 then bugmes();quit;end
 
//   fortran routine dres1 in dir. routines/default, without jocabian
 
yy5=dassl([y0,y0d],t0,t,atol,rtol,'dres1',info); //=yy4
 
if norm(dres1(yy5(1,:),yy5(2:3,:),yy5(4:5,:)),1)>1.d-5 then bugmes();quit;end
 
//   scilab functions,
 
yy6=dassl([y0,y0d],t0,t,atol,rtol,dres1,djac1,info);
 
if norm(dres1(yy6(1,:),yy6(2:3,:),yy6(4:5,:)),1)>1.d-5 then bugmes();quit;end
 
//   fortran routines
 
yy7=dassl([y0,y0d],t0,t,atol,rtol,'dres1','djac1',info); //==yy6
 
if norm(dres1(yy7(1,:),yy7(2:3,:),yy7(4:5,:)),1)>1.d-5 then bugmes();quit;end
 
 
//Copyright INRIA
 
 
//   Testing dassl/dasrt for E xdot - A x = 0  (index one)
 
//   x(0)=x0;   xdot(0)=xd0
 
Eps=1.e-9
 Eps  =
 
    1.000D-09  
 
rand('seed',0);
 
nx=5;
 
E=rand(nx,1)*rand(1,nx);A=rand(nx,nx);
 
//         Index 1
 
[Si,Pi,Di,o]=penlaur(E,A);pp=Si*E;[q,m]=fullrf(pp);
 
x0=q(:,1);x0d=pinv(E)*A*x0;
 
deff('[r,ires]=g(t,x,xdot)','r=E*xdot-A*x;ires=0');
 
t=[1,2,3];t0=0;info=list([],0,[],[],[],0,0);
 
x=dassl([x0,x0d],t0,t,g,info);x1=x(2:nx+1,:);
 
if norm(pp*x1-x1,1)  > Eps then bugmes();quit;end
 
//x(4) goes through 1 at  t=1.5409711:
 
t=1.5409711;ww=dassl([x0,x0d],t0,t,g,info);
 
if abs(ww(5)-1)  > 1000*Eps then bugmes();quit;end
 
//Find this fact with dasrt:
 
deff('[rt]=surface(t,y,yd)','rt=y(4)-1');nbsurf=1;
 
[yyy,nnn]=dasrt([x0,x0d],t0,t,g,nbsurf,surface,info);
 
//Defining Jacobian
 
deff('pd=j(t,y,ydot,cj)','pd=cj*D-A');
 
x=dassl([x0,x0d],t0,t,g,j,info);x2=x(2:nx+1,1);
 
if norm(x2-ww(2:nx+1,1),1)  > Eps then bugmes();quit;end
 
[yyy1,nnn]=dasrt([x0,x0d],t0,t,g,j,nbsurf,surface,info);
 
//Assuming x0d is not known:
 
x0d=ones(x0);info(7)=1;
 
x=dassl([x0,x0d],t0,t,g,info);
 
xn=dassl([x0,x0d],t0,t,g,j,info);
 
if norm(x-xn,1)  > Eps then bugmes();quit;end
 
 
//Copyright INRIA
 
Eps=1.e-2
 Eps  =
 
    0.01  
 
//     Example of use of ode function:
 
//     System to solve:
 
//     dy1/dt = -0.04*y1 + 1.e4*y2*y3
 
//     dy2/dt = 0.04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
 
//     dy3/dt = 3.e7*y2**2
 
// on the interval from t = 0.0 to t = 4.e10, with initial conditions
 
// y1 = 1.0, y2 = y3 = 0.  the problem is stiff.
 
//     Scilab definition of rhs function:
 
deff('[ydot]=f(t,y)',...
     'f1=-0.04*y(1)+1e4*y(2)*y(3),...
      f3=3e7*y(2)*y(2),...
      ydot=[f1;-f1-f3;f3]')
 
//     Definition of jacobian
 
deff('[jac]=j(t,y)',...
     'jac(1,1)=-0.04;jac(1,2)=1.e4*y(3);jac(1,3)=1.e4*y(2),...
      jac(3,1)=0;jac(3,2)=6.e7*y(2);jac(3,3)=0;...
      jac(2,1)=-jac(1,1);jac(2,2)=-jac(1,2)-jac(3,2);jac(2,3)=-jac(1,3);')
Warning :redefining function: j                       
 
 
//    y0= Initial state  t1= times at which sol. is computed
 
y0=[1;0;0];t0=0;t1=[0.4,4];nt=size(t1,'*');
 
//    Reference solution is:
 
yref=[0.9851721 0.9055180;0.0000339 0.0000224;0.0147940 0.0944596];
 
//
 
//  1. fortran program wfex and wjex dynamically called
 
files=G_make(["/tmp/wfex.o"],"wfex.dll");
 
iwfex=link(files,"wfex");
shared archive loaded
Link done
 
files=G_make(["/tmp/wjex.o"],"wjex.dll");
 
iwjex=link(files,"wjex");
shared archive loaded
Link done
 
//     jacobian is not given
 
y1=ode(y0,t0,t1,'wfex');
 
 
//  2. fortran called type given (stiff), no jacobian
 
y2=ode('stiff',y0,t0,t1,'wfex');
 
//  3. fortran called type (stiff) given
 
y3=ode('stiff',y0,t0,t1,'wfex','wjex');
 
//   hot restart
 
[z,w,iw]=ode('stiff',y0,0,0.4,'wfex','wjex');
 
z=ode('stiff',z,0.4,4,'wfex','wjex',w,iw);
 
if max(z-y3(:,2))  > Eps then bugmes();quit;end
 
 
[y1,w,iw]=ode(y0,t0,t1(1),'wfex');
 
y2=ode(y0,t1(1),t1(2:nt),'wfex',w,iw);
 
if max([y1 y2]-yref)  > Eps then bugmes();quit;end
 
 
[y1,w,iw]=ode(y0,t0,t1(1),'wfex','wjex');
 
y2=ode(y0,t1(1),t1(2:nt),'wfex','wjex',w,iw);
 
if max([y1 y2]-yref)  > Eps then bugmes();quit;end
 
 
// Unlink wfex and wjex
 
 
ulink(iwfex);
 
ulink(iwjex);
 
//   Setting tolerances
 
atol=[0.001,0.0001,0.001];rtol=atol;
 
//
 
//  4. type given , scilab lhs ,jacobian not passed
 
y4=ode('stiff',y0,t0,t1(1),atol,rtol,f);
 
if maxi(y4(:,1)-yref(:,1))  > Eps then bugmes();quit;end
 
 
//  5. type non given,  scilab rhs and jacobian functions
 
y5=ode(y0,t0,t1,f,j);
 
if max(y5-yref)  > Eps then bugmes();quit;end
 
//  6. type given (stiff),rhs and jacobian  by scilab
 
y6=ode('stiff',y0,t0,t1,0.00001,0.00001,f,j);
 
if max(y6-yref)  > Eps then bugmes();quit;end
 
 
//  7. matrix rhs, type given (adams),jacobian not passed
 
//
 
a=rand(3,3);ea=expm(a);
 
deff('[ydot]=f(t,y)','ydot=a*y')
Warning :redefining function: f                       
 
 
t1=1;y=ode('adams',eye(a),t0,t1,f);
 
if max(ea-y)  > Eps then bugmes();quit;end
 
 
 
 
 
//Copyright INRIA
 
Eps =1.e-5
 Eps  =
 
    0.00001  
 
//     Example of use of ode function:
 
//     System to solve:
 
//     dy1/dt = -0.04*y1 + 1.e4*y2*y3
 
//     dy2/dt = 0.04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
 
//     dy3/dt = 3.e7*y2**2
 
// on the interval from t = 0.0 to t = 4.e10, with initial conditions
 
// y1 = 1.0, y2 = y3 = 0.  the problem is stiff.
 
//     Scilab definition of rhs function:
 
deff('[ydot]=f(t,y)',...
     'f1=-0.04*y(1)+1e4*y(2)*y(3),...
      f3=3e7*y(2)*y(2),...
      ydot=[f1;-f1-f3;f3]')
Warning :redefining function: f                       
 
 
//     Definition of jacobian
 
deff('[jac]=j(t,y)',...
     'jac(1,1)=-0.04;jac(1,2)=1.e4*y(3);jac(1,3)=1.e4*y(2),...
      jac(3,1)=0;jac(3,2)=6.e7*y(2);jac(3,3)=0;...
      jac(2,1)=-jac(1,1);jac(2,2)=-jac(1,2)-jac(3,2);jac(2,3)=-jac(1,3);')
 
//    y0= Initial state  t1= times at which sol. is computed
 
y0=[1;0;0];t0=0;t1=[0.4,4];nt=size(t1,'*');
 
//    Reference solution is:
 
yref=[0.9851721 0.9055180;0.0000339 0.0000224;0.0147940 0.0944596];
 
//
 
//  Fortran program called in "routines/default" directory (Ex-ode.f)
 
//     jacobian is not given
 
y1=ode(y0,t0,t1,'fex');
 
if norm(yref-y1) > Eps then bugmes();quit;end
 
//     jacobian is given
 
y2=ode(y0,t0,t1,'fex','jex');
 
if norm(yref-y2) > Eps then bugmes();quit;end
 
 
 
 
 
 
//Copyright INRIA
 
 
// Dynamic call of fexab2 which reads A and B in Scilab's stack
 
deff('xd=f(t,x)','xd=A*x+B*sin(3*t)')
Warning :redefining function: f                       
 
 
A=rand(10,10)-4.5*eye();B=rand(10,1);
 
x=ode(ones(10,1),0,[1,2,3],f);
 
files=G_make(["/tmp/wfexab.o"],'wfexab.dll');
 
link(files,"wfexab");
shared archive loaded
Link done
 
if norm(x-ode(ones(10,1),0,[1,2,3],'wfexab')) > 10000*%eps then bugmes();quit;end
 
 
 
//Copyright INRIA
 
Eps=1.e-3
 Eps  =
 
    0.001  
 
deff('xcd=f(t,xc,xd,iflag)',...
    ['if iflag==0 then '
     '  xcd=fc(t,xc,e(t)-hd(t,xd));'
     'else '
     '  xcd=fd(xd,hc(t,xc));'
     'end']);
Warning :redefining function: f                       
 
 
 
A=[-10,2,3;4,-10,6;7,8,-10];B=[1;1;1];C=[1,1,1];
 
Ad=[1/2,1;0,1/20];Bd=[1;1];Cd=[1,1];
 
deff('st=e(t)','st=sin(3*t)')
 
deff('xdot=fc(t,x,u)','xdot=A*x+B*u')
 
deff('y=hc(t,x)','y=C*x')
 
deff('xp=fd(x,y)','xp=Ad*x + Bd*y')
 
deff('u=hd(t,x)','u=Cd*x')
 
 
h=0.1;t0=0;t=0:0.1:2;
 
x0c=[0;0;0];x0d=[0;0];nd=2;
 
xcd=odedc([x0c;x0d],nd,h,t0,t,f);
 
if norm(xcd-odedc([x0c;x0d],nd,h,t0,t,'fcd'),1)   > Eps then bugmes();quit;end
 
 
 
//(see default directory)
 
if norm(xcd-odedc([x0c;x0d],nd,h,t0,t,'fcd1'),1)  > Eps then bugmes();quit;end
 
 
deff('xcd=phis(t,xc,xd,iflag)',...
    ['if iflag==0 then '
     '  xcd=A*xc+B*xd;'
     'else '
     '  xcd=1-xd;'
     'end']);
 
 
t=0:0.1:30;
 
xcd=odedc([x0c;0],1,1,t0,t,phis);
 
//link('phis.o','phis')
 
xcd2=odedc([x0c;0],1,1,t0,t,'phis');
 
 
if norm(xcd-xcd2,1)  > Eps then bugmes();quit;end
 
deff('xd=ff(t,x)','xd=A*x+B*u')
 
 
u=1/2;xn=ode(x0c,t0,t,ff);
 
//plot2d([t',t',t',t'],[(xcd(1,:))',(xcd(2,:))',(xcd(3,:))',(xcd(4,:))'])
 
 
deff('xcd=phit(t,xc,xd,iflag)',...
    ['if iflag==0 then '
     '  xcd=[A*xc(1:3,:)+B*xc(4);xd];'
     'else '
     '  xcd=-xd;'
     'end']);
 
 
xcdt=odedc([x0c;1;1],1,1,t0,t,phit);
 
//link('phit.o','phit')
 
xcdt2=odedc([x0c;1;1],1,1,t0,t,'phit');
 
if norm(xcdt-xcdt2,1) > Eps then bugmes();quit;end
 
//plot2d([t',t',t',t'],[(xcdt(1,:))',(xcdt(2,:))',(xcdt(3,:))',(xcdt(4,:))'])
 
 
xcdt3=odedc('adams',[x0c;1;1],1,1,t0,t,'phit');
 
if norm(xcdt3-xcdt2,1)  > Eps then bugmes();quit;end
 
 
xcdt4=odedc('fix',[x0c;1;1],1,1,t0,t,'phit');
 
if norm(xcdt4-xcdt2,1) > Eps then bugmes();quit;end
 
 
xcdt5=odedc('stiff',[x0c;1;1],1,1,t0,t,'phit');
 
if norm(xcdt5-xcdt2,1) > Eps then bugmes();quit;end
 
 
 
//Copyright INRIA
 
Eps=1.e-5
 Eps  =
 
    0.00001  
 
//        %ODEOPTIONS
 
//
 
rand('seed',0);rand('normal');
 
nx=20;A=rand(nx,nx);A=A-4.5*eye();
 
deff('y=f(t,x)','y=A*x')
Warning :redefining function: f                       
 
 
deff('J=j(t,x)','J=A')
Warning :redefining function: j                       
 
 
x0=ones(nx,1);t0=0;t=[1,2,3,4,5];
 
nt=size(t,'*');
 
eAt=expm(A*t(nt));
 
 
//        Test itask=%ODEOPTIONS(1)
 
 
//itask=1  --->  usual call (t=vector of instants, solution at all t)
 
//========================
 
itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
 
jacflag=2;ml=-1;mu=-1;
 
%ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
 
xf=ode(x0,t0,t,f);   //lsoda
 
if norm(xf(:,nt)-eAt*x0)  > Eps then bugmes();quit;end
 
xfj=ode(x0,t0,t,f,j);   //lsoda with jacobian
 Warning: Jacobian external is given, but 
 not used!,  see %ODEOPTIONS(6)
 
 
if norm(xfj(:,nt)-eAt*x0)  > Eps then bugmes();quit;end
 
 
 
//itask=2;   --->  solution at mesh points  ---> t=tmax
 
//========================
 
%ODEOPTIONS(1)=2;tmax=t(5);
 
xft=ode(x0,t0,tmax,f);
 
[p,q]=size(xft);
 
xlast=xft(2:nx+1,q);
 
if xft(1,q)<tmax then bugmes();quit;end
 
if norm(xlast-expm(A*xft(1,q))*x0) > Eps  then bugmes();quit;end
 
 
//itask=3;   ---> solution at first mesh point beyond t=tmax
 
%ODEOPTIONS(1)=3;
 
x3=ode(x0,t0,tmax,f);
 
if norm(x3(2:nx+1)-xlast,1) > Eps then bugmes();quit;end
 
 
//itask=4; test with %ODEOPTIONS(2)=tcrit
 
%ODEOPTIONS(1)=4; //---> computation at all t and t>tcrit are not called
 
tcrit=2.5;%ODEOPTIONS(2)=tcrit;
 
chk=0;
 
deff('y=fcrit(t,x)',['if t<=tcrit then'
                      ' y=A*x;'
                      'else'
                      ' y=A*x;chk=resume(1);end'])
 
x42=ode(x0,t0,t,fcrit);
 Warning: integration up to tcrit
 
 
if chk==1 then bugmes();quit;end
 
[p,q]=size(x42);
 
if norm(x42(:,q)-ode(x0,t0,tcrit,f),1) > Eps then bugmes();quit;end
 
 
//itask=5; test with %ODEOPTIONS(2)=tcrit
 
%ODEOPTIONS(1)=5;  //---> computation at mesh points and t>tcrit are not called
 
%ODEOPTIONS(6)=2;  // Estimated jacobian
 
chk=0;
 
x52=ode(x0,t0,2.3,fcrit);
 
if chk==1 then bugmes();quit;end
 
[p,q]=size(x52);
 
if x52(1,q)>tcrit then bugmes();quit;end
 
 
//test of %ODEOPTIONS(3:5)=[h0,hmax,hmin]
 
%ODEOPTIONS(1)=1;
 
h0=0.0;hmax=0.1;hmin=0.0001;
 
%ODEOPTIONS(3:5)=[h0,hmax,hmin];
 
x35=ode(x0,t0,t,f);
 
if norm(x35-xf,1)  > Eps then bugmes();quit;end
 
 
//test of %ODEOPTIONS(6)=jacflag
 
%ODEOPTIONS(6)=1;//Jacobian given
 
%ODEOPTIONS(3:5)=[0 0 0];
 
x61=ode('st',x0,t0,t,f,j);   //with Jacobian
 
if norm (x61-xf,1) > Eps then bugmes();quit;end
 
 
 
%ODEOPTIONS(6)=0; // jacobian nor called nor estimated
 
x60=ode('st',x0,t0,t,f,j);   //Jacobian not used (warning)
 Warning: Jacobian external is given, but 
 not used!,  see %ODEOPTIONS(6)
 
 
x60=ode('st',x0,t0,t,f);    //Jacobian not used
 
if norm (x60-x61,1)  > Eps then bugmes();quit;end
 
 
 
%ODEOPTIONS(6)=1;//Jacobian estimated
 
x60=ode('st',x0,t0,t,f)  ;
 Warning: No Jacobian external given but 
 one is required by %ODEOPTIONS(6) value!
 
 
if norm (x60-x61,1) > Eps then bugmes();quit;end
 
 
//test of %ODEOPTIONS(6)=jacflag   (adams)
 
%ODEOPTIONS(6)=1;//with given Jacobian
 
x60=ode('ad',x0,t0,t,f,j) ;
 
if norm (x60-x61,1)  > Eps then bugmes();quit;end
 
 
 
%ODEOPTIONS(6)=0;// jacobian nor called nor estimated
 
x60=ode('ad',x0,t0,t,f,j);   //Jacobian not used (warning)
 Warning: Jacobian external is given, but 
 not used!,  see %ODEOPTIONS(6)
 
 
x60=ode('ad',x0,t0,t,f);    //Jacobian not used
 
if norm (x60-x61,1) > Eps then bugmes();quit;end
 
 
// test lsoda
 
%ODEOPTIONS(6)=2;// jacobian  estimated
 
x60=ode(x0,t0,t,f);
 
if norm (x60-x61,1) > Eps then bugmes();quit;end
 
 
%ODEOPTIONS(6)=1;//Jacobian estimated
 
x60=ode(x0,t0,t,f);
 Warning: No Jacobian external given but 
 one is required by %ODEOPTIONS(6) value!
 
 
if norm (x60-x61,1) > Eps then bugmes();quit;end
 
 
 
//   Banded Jacobian
 
itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
 
//provisional values as default
 
jacflag=2;ml=-1;mu=-1;
 
 
%ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
 
jacflag=2;%ODEOPTIONS(6)=jacflag;   //Banded Jacobian, given
 
nx=20;A=diag(-[1:nx]);x0=ones(nx,1);
 
t0=0;t=[1,2,3,4,5];
 
for k=1:nx-1, A(k,k+1)=1;end
 
for k=1:nx-2, A(k,k+2)=-2;end
 
for k=1:nx-1, A(k+1,k)=-1;end
 
for k=1:nx-2, A(k+2,k)=2;end
 
for k=1:nx-3, A(k+3,k)=-3;end
 
deff('xd=f(t,x)','xd=A*x')
Warning :redefining function: f                       
 
 
ml=3;mu=2;
 
%ODEOPTIONS(11:12)=[ml,mu];
 
for i=1:nx;
    for j=1:nx;
if A(i,j)<>0 then J(i-j+mu+1,j)=A(i,j);end
end;end;
Warning :redefining function: j                       
 
 
// J is a ml+mu+1 x ny matrix.
 
// Column 1 of J is made of mu zeros followed by df1/dx1, df2/dx1, df3/dx1,
 
// i.e. 1 + ml possibly nonzeros entries.
 
// Column 2 of J is made of mu-1 zeros followed by df1/dx2, df2/dx2,0...
 
// etc...
 
%ODEOPTIONS(6)=1;%ODEOPTIONS(11:12)=[-1,-1];
 
deff('jj=j1(t,x)','jj=A')
 
xnotband=ode('st',x0,t0,t,f,j1);
 
 
 
%ODEOPTIONS(6)=4;//banded jacobian external given
 
%ODEOPTIONS(11:12)=[3,2];
 
deff('jj=j(t,x)','jj=J')
 
xband=ode('st',x0,t0,t,f,j);
 
if norm (xnotband-xband,1)  > Eps then bugmes();quit;end
 
 
%ODEOPTIONS(6)=5;//banded jacobian evaluated
 
%ODEOPTIONS(11:12)=[3,2];
 
deff('jj=j(t,x)','jj=J')
 
xband=ode('st',x0,t0,t,f,j);
 Warning: Jacobian external is given, but 
 not used!,  see %ODEOPTIONS(6)
 
 
if norm (xnotband-xband,1) > Eps then bugmes();quit;end
 
 
 
//            Test of %ODEOPTIONS(7)
 
//%ODEOPTIONS(7)=mxstep  ---> maximum number od steps allowed
 
itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
 
//provisional values as default
 
jacflag=2;ml=-1;mu=-1;
 
%ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
 
%ODEOPTIONS(7)=10;
 
//ode(x0,t0,t,f);  // ---> Non convergence
 
 
//            Test of %ODEOPTIONS(8:9)
 
//%ODEOPTIONS(8:9)=[maxordn,maxords] ---> maximum order for nonstiff and stiff
 
itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
 
//provisional values as default
 
jacflag=2;ml=-1;mu=-1;
 
%ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
 
%ODEOPTIONS(8:9)=[0,0]; //---> default values
 
wref=ode(x0,t0,t,f); //just for computing reference
 
 
%ODEOPTIONS(8:9)=[4,3];
 
ww=ode(x0,t0,t,f);
 
if norm (wref-ww,1) > Eps then bugmes();quit;end
 
 
%ODEOPTIONS(8:9)=[12,5];
 
if norm (wref-ode(x0,t0,t,f),1) > Eps then bugmes();quit;end
 
 
//using stiff method
 
 
%ODEOPTIONS(9)=0;
 
wref=ode('st',x0,t0,t,f);
 
%ODEOPTIONS(9)=5;
 
if norm (wref-ode('st',x0,t0,t,f),1)  > Eps then bugmes();quit;end
 
%ODEOPTIONS(9)=4;
 
if norm (wref-ode('st',x0,t0,t,f),1) > Eps then bugmes();quit;end
 
 
 
//using nonstiff method
 
 
%ODEOPTIONS(8)=0;
 
wref=ode('ad',x0,t0,t,f);
 
%ODEOPTIONS(8)=12;
 
if norm (wref-ode('ad',x0,t0,t,f),1)  > Eps then bugmes();quit;end
 
%ODEOPTIONS(8)=5;
 
if norm (wref-ode('ad',x0,t0,t,f),1) > Eps then bugmes();quit;end
 
 
//mixed
 
%ODEOPTIONS(8:9)=[5,12];
 
wref=ode(x0,t0,t,f);
 
%ODEOPTIONS(8:9)=[4,10];
 
if norm (ode(x0,t0,t,f)-wref,1) > Eps then bugmes();quit;end
 
 
 
A=diag([-10,-0.01,-1]);
 
deff('uu=u(t)','uu=sin(t)');
 
B=rand(3,1);
 
deff('y=f(t,x)','y=A*x+B*u(t)')
Warning :redefining function: f                       
 
 
%ODEOPTIONS(1)=2;
 
yy1=ode('stiff',[1;1;1],0,1,f);
 
yy2=ode('stiff',[1;1;1],0,2,f);
 
%ODEOPTIONS(1)=3;
 
yy1=ode('stiff',[1;1;1],0,1,f);
 
yy2=ode('stiff',[1;1;1],0,2,f);
 
 
clear %ODEOPTIONS;
 
rand('seed',0);rand('normal');
 
nx=20;A=rand(nx,nx);A=A-4.5*eye();
 
deff('y=f(t,x)','y=A*x')
Warning :redefining function: f                       
 
 
deff('J=j(t,x)','J=A')
Warning :redefining function: j                       
 
 
//%ODEOPTIONS(1)=1;
 
y2=ode('stiff',ones(nx,1),0,2,f,j);
 
[y1,w,iw]=ode('stiff',ones(nx,1),0,1,f,j);
 
y2p=ode('stiff',y1,1,2,f,j,w,iw);
 
y12=ode('stiff',ones(nx,1),0,[1,2],f,j);
 
if norm (y12(:,2)-y2p) > Eps then bugmes();quit;end
 
yaf=ode('adams',ones(nx,1),0,2,f,j);
 
yaj=ode('adams',ones(nx,1),0,2,f,j);
 
ysf=ode('stiff',ones(nx,1),0,2,f,j);
 
ysj=ode('stiff',ones(nx,1),0,2,f,j);
 
 
 
 
