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

subs={b[n_]-> (Subscript[b, n]),d[n_]-> (Subscript[d, n]),k[n_]-> (Subscript[k, n]),q[n_]-> (Subscript[q, n]),x[n_]->(Subscript[x, n]),x'[n_]->Subsuperscript[x,n,′],b'[n_]->Subsuperscript[b,n,′],k'[n_]->Subsuperscript[k,n,′],q'[n_]->Subsuperscript[q,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

P[D[x'[1],x'[0]]]

(* 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]+k[n]dx[n][i]/;n≥1 && i≤2

dx'[n_][3]:=dx'[n-1][3]+q[n]+k[n]dx[n][3]/; n≥1

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

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

P[Simplify[dx'[2][3]]]

(* The determinant for closed-orbit evolution *)

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

P[Simplify[cD[2]]]

(* Horizontal phase advance value *)

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

P[Simplify[twocosphiX[2]]]

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])

P[Simplify[dxc[2]]]

P[Simplify[dxcp[2]]]

(* 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]+k'[n](dx[n][1])^2+k[n]ddx[n][1]/;n≥1

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

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

ddx'[n_][4]:=ddx'[n-1][4]+q'[n]dx[n][1]+k'[n]dx[n][1]dx[n][3]+k[n]ddx[n][4]/;n≥1

ddx'[n_][5]:=ddx'[n-1][5]+q'[n]dx[n][2]+k'[n]dx[n][2]dx[n][3]+k[n]ddx[n][5]/;n≥1

P[Array[ddx'[1],5]]

P[Array[ddx[2],5]]

P[Simplify[Array[ddx'[2],5]]]

(* 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[Simplify[dtwocosphiX[2]]]

(* 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]-k[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_]:=q'[n]+k'[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]-k[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]]]