function Xdot=f(t,X)
P=X(1)
A=X(2)
Z=X(3)
D=X(4)
//------------------------------------------------------
Xmat1=[
0, A*P/(1+P), 0, 0;
0, 0, 0.01*Z*0.1, 0.01*Z*0.9 + 0.01*A;
0, 0, 0, 0.1*Z;
0.01*D, 0, 0, 0
]
dX1=[
-sum( Xmat1(1,:) ) + sum( Xmat1(:,1) );
-sum( Xmat1(2,:) ) + sum( Xmat1(:,2) );
-sum( Xmat1(3,:) ) + sum( Xmat1(:,3) );
-sum( Xmat1(4,:) ) + sum( Xmat1(:,4) );
]
//------------------------------------------------------
P=X(5)
A=X(6)
Z=X(7)
D=X(8)
//------------------------------------------------------
Xmat2=[
0, A*P/(1+P), 0, 0;
0, 0, 0.01*Z*0.1, 0.01*Z*0.9 + 0.01*A;
0, 0, 0, 0.1*Z;
0.01*D, 0, 0, 0
]
dX2=[
-sum( Xmat2(1,:) ) + sum( Xmat2(:,1) );
-sum( Xmat2(2,:) ) + sum( Xmat2(:,2) );
-sum( Xmat2(3,:) ) + sum( Xmat2(:,3) );
-sum( Xmat2(4,:) ) + sum( Xmat2(:,4) );
]
Xdot=[dX1 * 0.9 + dX2 * 0.1;
dX1 * 0.1 + dX2 * 0.9]
endfunction
// 初期値
X1s=[2; 1; 0.01; 0]
X2s=[2; 3; 0.01; 1]
Xs=[X1s; X2s]
//時間
t0=0; t1=50; dt=0.01;
t= t0 : dt : t1;
//微分方程式solver
x=ode(Xs,t0,t,f)
// グラフの表示
clf()
scf(0)
plot2d(t,x(1,:),2)
plot2d(t,x(2,:),3)
plot2d(t,x(3,:),4)
plot2d(t,x(4,:),5)
scf(1)
plot2d(t,x(5,:),6)
plot2d(t,x(6,:),7)
plot2d(t,x(7,:),8)
plot2d(t,x(8,:),9)
