g(d):=δ/2 * d;

Li(si,sj,d) := a_4*exp(a_1*si + a_2 * sj + a_3*d)$
Si_tt(si,sj,d,xi) :=  xi + si*(1-Li(si,sj,d))$
D_tt(si,sj,dt,xi,xj):= (1-δ)*dt + g(dt) + Γ*(xi + xj) + γ *(Li(si,sj,dt) + Li(sj,si,dt))$

B: matrix(
 [1,0,Γ], 
 [0,1,Γ]
)$

A(si,sj,dt,xi,xj) := transpose(matrix(
    [diff(Si_tt(si,sj,dt,xi),si),  diff(Si_tt(si,sj,dt,xi),sj), diff(Si_tt(si,sj,dt,xi),dt) ],
    [diff(Si_tt(sj,si,dt,xi),si),  diff(Si_tt(sj,si,dt,xi),sj), diff(Si_tt(sj,si,dt,xi),dt) ],
    [diff(D_tt(si,sj,dt,xi,xj),si),  diff(D_tt(si,sj,dt,xi,xj),sj),  diff(D_tt(si,sj,dt,xi,xj),dt)]
))$

BAt : B . invert(A(s1,s2,d,x1,x2)) ;

The next step is to create an iterated version.

help()subst)
;