//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']);
 
 
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);
 
 
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);
 
 
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