u[A] := delta*k*(P[1]*cosh(k*y) + P[2]*k*y*sinh(k*y))*sin(k*x);
v[A] := (1+nu)*y*eta/(1-nu)+delta*k*(((3-4*nu)*P[2]-P[1])*sinh(k*y) - P[2]*k*y*cosh(k*y))*cos(k*x);
epsilon[1] := diff(u[A],x);
epsilon[2] := diff(v[A],y);
epsilon[3] := 0.5*(diff(u[A],y) + diff(v[A],x));
sigma[1] := E[A]/((1+nu)*(1-2*nu))*((1-nu)*(epsilon[1]-eta*(1+nu)) + nu*(epsilon[2]-eta*(1+nu)));
sigma[2] := E[A]/((1+nu)*(1-2*nu))*((1-nu)*(epsilon[2]-eta*(1+nu)) + nu*(epsilon[1]-eta*(1+nu)));
sigma[3] := E[A]/(1+nu) *epsilon[3];
w[1] := simplify((sigma[1]*(epsilon[1]-eta) + sigma[2]*(epsilon[2]-eta) + 2*sigma[3]*epsilon[3])/2);
u[B] := delta*k*(M[1]-M[2]*k*y)*exp(-k*y)*sin(k*x);
v[B] := h*(1+nu)*eta/(2*(1-nu)) + delta*k*((M[1] -(3-4*nu)*M[2])-M[2]*k*y)*exp(-k*y)*cos(k*x);
epsilon[4] := diff(u[B],x);
epsilon[5] := diff(v[B],y);
epsilon[6] := 0.5*(diff(u[B],y) + diff(v[B],x));
sigma[4] := E[B]/((1+nu)*(1-2*nu))*((1-nu)*(epsilon[4]) + nu*(epsilon[5]));
sigma[5] := E[B]/((1+nu)*(1-2*nu))*((1-nu)*(epsilon[5]) + nu*(epsilon[4]));
sigma[6] := E[B]/(1+nu) *epsilon[6];
w[2] := simplify((sigma[4]*epsilon[4] + sigma[5]*epsilon[5] + 2*sigma[6]*epsilon[6])/2);
w[3] := simplify(sigma[4]*(epsilon[4]-epsilon[1]) + sigma[5]*(epsilon[5]-epsilon[2]) + sigma[6]*(epsilon[6]-epsilon[3]));
w[4] := simplify(-Einterfacial*delta*k*k*cos(k*x) + w[1] - w[2] + w[3]);
w[5] := diff(w[4],x);
w[6] := diff(w[5],x);
with(linalg):
U := matrix(4,4,[-cosh(omega), -omega*sinh(omega), (cosh(omega)-sinh(omega)),(sinh(omega)-cosh(omega))*omega,sinh(omega), omega*cosh(omega)-5*sinh(omega)/3,(cosh(omega)-sinh(omega)),(sinh(omega)-cosh(omega))*(5+3*omega)/3,cosh(omega),(3*omega*sinh(omega)-4*cosh(omega))/3,(1-alpha)*(sinh(omega)-cosh(omega))/(1+alpha),(1-alpha)*(cosh(omega)-sinh(omega))*(4+3*omega)/(3*(1+alpha)),sinh(omega),(3*omega*cosh(omega)-sinh(omega))/3,(1-alpha)*(cosh(omega)-sinh(omega))/(1+alpha),(1-alpha)*(sinh(omega)-cosh(omega))*(1+3*omega)/(3*(1+alpha))]);
V := matrix(4,1,[0,2.*eta/k,0,2*eta/k]);
Uinv := inverse(U);
PM := evalm(Uinv &* V);
P[1] := simplify(PM[1,1]);
P[2] := simplify(PM[2,1]);
M[1] := simplify(PM[3,1]);
M[2] := simplify(PM[4,1]);
w[7] := taylor(w[6],delta=0,2);
w[8] := subs(nu=1./3.,y=omega/k,E[B] = E[A]*(1-alpha)/(1+alpha),w[7]);
w[9] := simplify(w[8]);