#
#To run this program, load it into a file, for example `t'
#
#Type `xmaple'
#
#When the screen appears, type `read t;'
#
#Make sure you put the final `;'
#
#The program will then run in the screen. 
#
#Notice that anything following the symbol # on the same line is 
#disregarded by Maple.
#
#
phi:=C1*x^2+C2*x*y+C3*y^2+C4*x^3+C5*x^2*y+C6*x*y^2+C7*y^3+C8*x^4+C9*x^3*y+C10*x^2*y^2+C11*x*y^3+C12*y^4;
#
#
sxx1:= diff(phi,y,y); 
syy1:= diff(phi,x,x); 
sxy1:= -diff(phi,x,y); 
#
#The command `unapply(...,x,y) configures the stresses as functions of x,y,
#so that we can find the value at (e.g.) y=b.
#
sxx2:=unapply(sxx1,x,y);
syy2:=unapply(syy1,x,y);
sxy2:=unapply(sxy1,x,y);
#
#We now find the tractions on y=b as
#
t1:=syy2(x,b);
t2:=sxy2(x,b);
#
#and on y=-b
#
t3:=syy2(x,-b);
t4:=sxy2(x,-b);
#
#On x=0, we have
#
t5:=sxx2(0,y);
t6:=sxy2(0,y);
#
#The stress function is order 4, so the stresses are order 2 in x and y.
#The tractions on y=+b or -b might therefore be polynomials in x of order 2.
#
#We calculate the coefficients of each power of x in these expressions as
#
s1:=coeff(t1,x,2);
s2:=coeff(t1,x,1);
s3:=coeff(t1,x,0);
s4:=coeff(t2,x,2);
s5:=coeff(t2,x,1);
s6:=coeff(t2,x,0);
s7:=coeff(t3,x,2);
s8:=coeff(t3,x,1);
s9:=coeff(t3,x,0);
s10:=coeff(t4,x,2);
s11:=coeff(t4,x,1);
s12:=coeff(t4,x,0);
#
#The biharmonic equation is 4th order, so applying it to a 4th order polynomial
#generates a constant.
#
biharm:=diff(phi,x$4)+diff(phi,y$4)+2*diff(phi,x,x,y,y);
#
#and this must be a constant which has to be zero.
#
#We also calculate the three force resultants on x=0 as
#
Fx:=int(t5, y=-b..b);
Fy:=int(t6, y=-b..b);
M:=int(t5*y, y=-b..b);
#
#We now solve for the constants so as to satisfy (i) the strong boundary
#conditions, (ii) the biharmonic equation and (iii) the weak boundary
#conditions.
#
solution:=solve({s1=0,s2=0,s3=0,s4=0,s5=0,s6=0,s7=0,s8=0,s9=0,s10=0,s11=0,s12=0,biharm=0,Fx=0,M=0,Fy=F},{C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12});
#
#Notice that there are more equations than there are constants. Some of the 
#equations are not linearly independent. However, Maple can handle this if
#there is a solution.
#
#Substitute the solution into the original stress function.
#
phi:=subs(solution,phi);
#
#and calculate the final stresses as
#
sxx3:=diff(phi,y,y); 
syy3:=diff(phi,x,x); 
sxy3:=-diff(phi,x,y); 
#
#and for fun, make a contour plot of (say) sxx3
#
sxx4:=unapply(sxx3,x,y,F,b);
#
#Notice that I had to make it a function of x,y and F,b so that I could put
#numbers in for the latter variables.
#
with(plots):
contourplot(sxx4(x,y,1,1),x=0..10,y=-1..1,filled=true,coloring=[red,blue]);
#
#Notice that you can make a scale picture of the body from this plot by
#clicking inside the picture with the first mouse button and then clicking
#on 1:1 on the pop up tool bar. Otherwise the scale used for x and y will
#be different.