Timoshenko beam

import sympy as sp
w, phi = sp.symbols('w phi', cls=sp.Function)
C1, C2, C3, C4 = sp.symbols('C1 C2 C3 C4')
x, EI, GA, L, q = sp.symbols('x EI GA L q')

diffeq1 = sp.Eq(EI*sp.diff(phi(x),x,2)-GA*(sp.diff(w(x),x,1)+phi(x)),0)
diffeq2 = sp.Eq(GA*(sp.diff(w(x),x,2)+sp.diff(phi(x),x,1)),-q)
display(diffeq1)
display(diffeq2)

w,phi = sp.dsolve([diffeq1, diffeq2],[w(x),phi(x)])
w = w.rhs
phi = phi.rhs

\[\displaystyle EI \frac{d^{2}}{d x^{2}} \phi{\left(x \right)} - GA \left(\phi{\left(x \right)} + \frac{d}{d x} w{\left(x \right)}\right) = 0\]

\[\displaystyle GA \left(\frac{d}{d x} \phi{\left(x \right)} + \frac{d^{2}}{d x^{2}} w{\left(x \right)}\right) = - q\]

gamma_shear = phi + sp.diff(w, x)
kappa = sp.diff(phi, x)
M = EI * kappa
V = GA*gamma_shear

eq1  = sp.Eq(w.subs(x , 0) , 0)
eq2  = sp.Eq(M.subs(x , 0) , 0)
eq3  = sp.Eq(w.subs(x , L) , 0)
eq4  = sp.Eq(M.subs(x , L) , 0)

sol = sp.solve((eq1,eq2,eq3,eq4) ,
               (C1 ,C2 ,C3 ,C4))
w_sol = w.subs(sol)
M_sol = M.subs(sol)

display(sp.simplify(w_sol.subs(x,L/2)))

\[\displaystyle \frac{L^{2} q}{8 GA} + \frac{5 L^{4} q}{384 EI}\]