# phs8x8.ode
# relative phases 
# 8x8 phase model
# row 1
x1_1p=w+h(x1_2)+h(x2_1)
x1_1'=-x1_1
x1_[2..7]'=w-x1_1p+h(x1_[j-1]-x1_[j])+h(x1_[j+1]-x1_[j])+h(x2_[j]-x1_[j])
x1_8'=w-x1_1p+h(x1_7-x1_8)+h(x2_8-x1_8)
# row 2
x2_1'=w-x1_1p+h(x1_1-x2_1)+h(x3_1-x2_1)+h(x2_2-x2_1)
x2_[2..7]'=w-x1_1p+h(x2_[j-1]-x2_[j])+h(x2_[j+1]-x2_[j])+h(x1_[j]-x2_[j])+h(x3_[j]-x2_[j])
x2_8'=w-x1_1p+h(x2_7-x2_8)+h(x1_8-x2_8)+h(x3_8-x2_8)
# row 3
x3_1'=w-x1_1p+h(x2_1-x3_1)+h(x4_1-x3_1)+h(x3_2-x3_1)
x3_[2..7]'=w-x1_1p+h(x3_[j-1]-x3_[j])+h(x3_[j+1]-x3_[j])+h(x2_[j]-x3_[j])+h(x4_[j]-x3_[j])
x3_8'=w-x1_1p+h(x3_7-x3_8)+h(x2_8-x3_8)+h(x4_8-x3_8)
# row 4
x4_1'=w-x1_1p+h(x3_1-x4_1)+h(x5_1-x4_1)+h(x4_2-x4_1)
x4_[2..7]'=w-x1_1p+h(x4_[j-1]-x4_[j])+h(x4_[j+1]-x4_[j])+h(x3_[j]-x4_[j])+h(x5_[j]-x4_[j])
x4_8'=w-x1_1p+h(x4_7-x4_8)+h(x3_8-x4_8)+h(x5_8-x4_8)
# row 5
x5_1'=w-x1_1p+h(x4_1-x5_1)+h(x6_1-x5_1)+h(x5_2-x5_1)
x5_[2..7]'=w-x1_1p+h(x5_[j-1]-x5_[j])+h(x5_[j+1]-x5_[j])+h(x4_[j]-x5_[j])+h(x6_[j]-x5_[j])
x5_8'=w-x1_1p+h(x5_7-x5_8)+h(x4_8-x5_8)+h(x6_8-x5_8)
# row 6
x6_1'=w-x1_1p+h(x5_1-x6_1)+h(x7_1-x6_1)+h(x6_2-x6_1)
x6_[2..7]'=w-x1_1p+h(x6_[j-1]-x6_[j])+h(x6_[j+1]-x6_[j])+h(x5_[j]-x6_[j])+h(x7_[j]-x6_[j])
x6_8'=w-x1_1p+h(x6_7-x6_8)+h(x5_8-x6_8)+h(x7_8-x6_8)
# row 7
x7_1'=w-x1_1p+h(x6_1-x7_1)+h(x8_1-x7_1)+h(x7_2-x7_1)
x7_[2..7]'=w-x1_1p+h(x7_[j-1]-x7_[j])+h(x7_[j+1]-x7_[j])+h(x6_[j]-x7_[j])+h(x8_[j]-x7_[j])
x7_8'=w-x1_1p+h(x7_7-x7_8)+h(x6_8-x7_8)+h(x8_8-x7_8)
# row 8
x8_1'=w-x1_1p+h(x8_2-x8_1)+h(x7_1-x8_1)
x8_[2..7]'=w-x1_1p+h(x8_[j-1]-x8_[j])+h(x8_[j+1]-x8_[j])+h(x7_[j]-x8_[j])
x8_8'=w-x1_1p+h(x8_7-x8_8)+h(x7_8-x8_8)
par a0=0,a1=0,a2=0,b1=1,b2=0,w=1
h(u)=a0+a1*cos(u)+a2*cos(2*u)+b1*sin(u)+b2*sin(2*u)
@ total=100,meth=euler
done