clear all
syms b rho h x y v PI h EE nu q U
n_terms = 3;
rho = 2;
a = sym('A',[n_terms,n_terms]);
v = 0;
for i=1:n_terms
    for j=1:n_terms
       v=v+a(i,j)*sin((2*i-1)*pi*x/b)*sin((2*j-1)*pi*y/(rho*b));
    end
end
U = simplify(EE*h^3/(24*(1-nu))*( (diff(v,x,2)+diff(v,y,2))^2 ...
    - 2*(1-nu)*( diff(v,x,2)*diff(v,y,2) - diff(diff(v,x),y)^2 ) ));
PI = int(int(U - q*v,x,[0,b]),y,[0,rho*b]);
count = 1;
for i = 1:n_terms
    for j = 1:n_terms
        eq(count) = diff(PI,a(i,j))==0;
        unknowns(count) = a(i,j);
        count = count + 1;
    end
end
asol = struct2cell(solve(eq,unknowns));

 for i=1:n_terms*n_terms
     v = subs(v,unknowns(i),asol{i});
 end
 simplify(v)
 fsurf(subs(v,[EE,nu,h,b,q],[1,0,1,1,1]),[0 1 0 rho])
 daspect([1 1 1])
 xlabel('x/b');
 zlabel('vEh^3/(1-\nu)b^4q');
 ylabel('y/\rho b');