(* Reformatting rules to make output nicer if you use /.subs *)

subs={b[n_]-> (Subscript[b, n]),d[n_]-> (Subscript[d, n]),x[n_]->(Subscript[x, n]),x'[n_]->Subsuperscript[x,n,′],b'[n_]->Subsuperscript[b,n,′],b''[n_]->Subsuperscript[b,n,′′]}

P[x_]:=TraditionalForm[x/.subs]

(* The basic dynamics, but don't want to expand this fully *)

x[n_]:=x[n-1]+d[n]x'[n-1]/;n≥1

x'[n_]:=x'[n-1]+p^-1b[n][x[n]]/;n≥1

D[x'[1],x'[0]] /.subs

(* First derivatives in my more compact notation *)

dx[n_][i_]:=dx[n-1][i]+d[n]dx'[n-1][i]/;n≥1

dx'[n_][i_]:=dx'[n-1][i]+p^-1b'[n]dx[n][i]/;n≥1 && i≤2

dx'[n_][3]:=dx'[n-1][3]-p^-2b[n]+p^-1b'[n]dx[n][3]/; n≥1

dx[0][i_]:=If[i==1,1,0]

dx'[0][i_]:=If[i==2,1,0]

Simplify[dx'[2][3]]/.subs

(* The determinant for closed-orbit evolution *)

cD[N_]:=(dx[N][1]-1)(dx'[N][2]-1)-dx[N][2]dx'[N][1]

Simplify[cD[2]]/.subs

(* Horizontal phase advance value *)

twocosphiX[N_]:=dx[N][1]+dx'[N][2]

Simplify[twocosphiX[2]]/.subs

Simplify[twocosphiX[2]+cD[2]]

Simplify[twocosphiX[3]+cD[3]]

Simplify[twocosphiX[N]+cD[N]]

(* Closed orbit evolution *)

dxc[N_]:=(1/cD[N])((1-dx'[N][2])dx[N][3]+dx[N][2]dx'[N][3])

dxcp[N_]:=(1/cD[N])(dx'[N][1]dx[N][3]+(1-dx[N][1])dx'[N][3])

Simplify[dxc[2]]/.subs

Simplify[dxcp[2]]/.subs

(* Second derivatives of beamline function:
1=xx 2=xx' 3=x'x' 4=xp 5=x'p *)

ddx[0][i_]:=0

ddx'[0][i_]:=0

ddx[n_][i_]:=ddx[n-1][i]+d[n]ddx'[n-1][i]/;n≥1

ddx'[n_][1]:=ddx'[n-1][1]+p^-1b''[n](dx[n][1])^2+p^-1b'[n]ddx[n][1]/;n≥1

ddx'[n_][2]:=ddx'[n-1][2]+p^-1b''[n](dx[n][1])(dx[n][2])+p^-1b'[n]ddx[n][2]/;n≥1

ddx'[n_][3]:=ddx'[n-1][3]+p^-1b''[n](dx[n][2])^2+p^-1b'[n]ddx[n][3]/;n≥1

ddx'[n_][4]:=ddx'[n-1][4]-p^-2b'[n]dx[n][1]+p^-1b''[n]dx[n][1]dx[n][3]+p^-1b'[n]ddx[n][4]/;n≥1

ddx'[n_][5]:=ddx'[n-1][5]-p^-2b'[n]dx[n][2]+p^-1b''[n]dx[n][2]dx[n][3]+p^-1b'[n]ddx[n][5]/;n≥1

ddx'[1][3]/.subs

(* Horizontal phase advance evolution *)

dtwocosphiX[N_]:=(ddx[N][1]+ddx'[N][2])dxc[N]+(ddx[N][2]+ddx'[N][3])dxcp[N]+ddx[N][4]+ddx'[N][5]

P[f=Simplify[dtwocosphiX[2]]]

P[f=Simplify[p^5*dtwocosphiX[2]*cD[2]]]

P[m=Simplify[CoefficientList[f,{b''[1],b''[2]}]]]

P[m[[1]][[1]]]

P[m[[1]][[2]]]

P[m[[2]][[1]]]

m[[2]][[2]]

sf:=(m[[1]][[1]]+m[[2]][[1]]b''[1]+m[[1]][[2]]b''[2])

P[sf]

Simplify[f-sf]

(* First derivatives of the Y plane beamline function *)

dy[n_][i_]:=dy[n-1][i]+d[n]dy'[n-1][i]/;n≥1

dy'[n_][i_]:=dy'[n-1][i]-p^-1b'[n]dy[n][i]/;n≥1

dy[0][i_]:=If[i==1,1,0]

dy'[0][i_]:=If[i==2,1,0]

P[Simplify[dy'[2][1]]]

(* Vertical phase advance value *)

twocosphiY[N_]:= dy[N][1]+dy'[N][2]

P[Simplify[twocosphiY[2]]]

(* Second Y plane derivatives: i=1 dy0dp, i=2 dy'0dp *)

ddy[0,N_][i_]:=0

ddy'[0,N_][i_]:=0

ddy[n_,N_][i_]:=ddy[n-1,N][i]+d[n]ddy'[n-1,N][i]/;n≥1

tdx[n_,N_]:=dx[n][1]dxc[N]+dx[n][2]dxcp[N]+dx[n][3]

P[Simplify[tdx[0,2]-tdx[2,2]]]

tdk[n_,N_]:=-p^-2b'[n]+p^-1b''[n]tdx[n,N]

P[Simplify[Array[tdk,{2,1},{1,2}]]]

ddy'[n_,N_][i_]:=ddy'[n-1,N][i]-tdk[n,N]dy[n][1]-p^-1b'[n]ddy[n,N][i]/;n≥1

P[Simplify[Array[ddy'[1,2],2]]]

dtwocosphiY[N_]:=ddy[N,N][1]+ddy'[N,N][2]

P[Simplify[dtwocosphiY[2]]]

P[f=Simplify[dtwocosphiY[2]*cD[2]*p^5]]

P[m=Simplify[CoefficientList[f,{b''[1],b''[2]}]]]

sf:=(m[[1]][[1]]+m[[2]][[1]]b''[1]+m[[1]][[2]]b''[2])

P[sf]

P[Simplify[sf-f]]