restart;
eq[1] := (D[2](U))(x, t+(1/2)*dt) = (D[1, 1](U))(x, t+(1/2)*dt);

eq[2] := (D(unapply(eq[1], x)))(x);

eq[3] := (D(unapply(eq[1], t)))(t);

eq[4] := (D(unapply(eq[2], x)))(x);

eq[5] := (D(unapply(eq[3], t)))(t);

eq[6] := (D(unapply(eq[3], x)))(x);

eq[7] := (D(unapply(eq[4], t)))(t);

eq[8] := (D(unapply(eq[4], x)))(x);

eq[9] := (D(unapply(eq[8], x)))(x);

eq[10] := (D(unapply(eq[9], t)))(t);

eq[11] := (D(unapply(eq[9], x)))(x);

eq[12] := (D(unapply(eq[11], x)))(x);

A1 := series(U(x, t+dt), dt = (1/2)*h, 4);

A2:=convert(subs(h=dt,A1),polynom);

A3 := subs({eq[1], eq[3], eq[5]}, A2);

A4 := subs({eq[4], eq[7]}, A3);

A5 := subs(eq[9], A4);

A6 := series(U(x, t), t = h, 4);

A7 := convert(subs(h = t+(1/2)*dt, A6), polynom);

A8 := subs(eq[9], subs({eq[4], eq[7]}, subs({eq[1], eq[3], eq[5]}, A7)));

Left := simplify((A5-A8)/dt);

B1 := series(U(x-dx, t+dt)-2*U(x, t+dt)+U(x+dx, t+dt), dx = 0);

B2 := convert(B1, polynom);

B3 := convert(subs(h = dt, series(B2, dt = (1/2)*h, 4)), polynom);

B4 := subs(eq[12], subs({eq[10], eq[12]}, subs(eq[9], subs({eq[7], eq[9], eq[4], eq[8]}, B3))));

C1 := convert(series(U(x-dx, t)-2*U(x, t)+U(x+dx, t), dx = 0), polynom);

C2 := convert(subs(h = t+(1/2)*dt, series(C1, t = h, 4)), polynom);

C3 := subs(eq[12], subs({eq[10], eq[12]}, subs(eq[9], subs({eq[7], eq[9], eq[4], eq[8]}, C2))));

Right := simplify((theta*B4+(1-theta)*C3)/dx^2);

T := collect(simplify(Left-Right), dt);

# Crank Nicholson Method
subs(theta = 1/2, T);

# 4th order in space
simplify(subs(theta = 1/2-dx^2/(12*dt), T));