kill(all);

elastic isotropic material plane stress relation between 
in-plane stress/strain components

D:
E/(1-nu^2)*
matrix(
[ 1,nu,       0],
[nu, 1,       0],
[ 0, 0,(1-nu)/2]
);

laminate generalized stress/strain component relation,
expressed as the three A,B and C stiffness matrix terms
here the laminate consists in a uniform material plate.

membrane stiffness term

A: integrate(D,z,-h/2,h/2);

membrane/bending coupling stiffness term

B: integrate(D*z,z,-h/2,h/2);

bending stiffness term

C: integrate(D*z^2,z,-h/2,h/2);

generalized stress/strain components for the laminate

e:matrix([ex],[ey],[gxy])$
k:matrix([kx],[ky],[kxy]);
q:matrix([qx],[qy],[qxy])$
m:matrix([mx],[my],[mxy])$

a usage example for the "matrix()" Maxima instruction

matrix( [1,2,3], [4,5,6]);

generalized stress/strain relations associated to moment resultants,
expressed in terms of a residual vector

m - B.e - C.k;

the terms of the residual vector are collected into 
a list of residuals for compatibility with the linsolve instruction

eqns:list_matrix_entries(%);

reference bending moment flow component, 
based on the beam theory, i.e. 
F*l bending moment uniformly distributed over a b width extension 

mref:F*l/b;

reference curvature, from the beam theory

kref:12*F*l/E/b/h^3;

Case A: zero "my" moment flow component, free ky curvature

null_my:[mx=mref,my=0,mxy=0];

null_my:append(
    null_my,
    linsolve(ev(eqns,null_my),[kx,ky,kxy])
);

Case A results

fullratsimp( 
    ev( 
[mx,my,mxy]/mref   , null_my));%,nu=0.3;

fullratsimp( 
    ev( 
[kx,ky,kxy]/kref   , null_my));%,nu=0.3;

Case B: zero ky curvature, "my" moment flow component to be determined.

null_ky:[mx=mref,ky=0,mxy=0];

null_ky:append(
    null_ky,
    linsolve(ev(eqns,null_ky),[kx,my,kxy])
);

Case B results

fullratsimp( 
    ev( 
[mx,my,mxy]/mref, null_ky));%,nu=0.3;

fullratsimp( 
    ev( 
[kx,ky,kxy]/kref, null_ky));%,nu=0.3;