--> //initializations

--> c=20 //KN/m^2 (cohesion)
 c  = 

   20.


--> x=12 // (angle of shearin resistance)
 x  = 

   12.


--> Y=18 //KN/m^3 (unit weight)
 Y  = 

   18.


--> z1=0 //m (at top)
 z1  = 

   0.


--> z2=4 //m (at bottom)
 z2  = 

   4.


--> //calculations

--> Ka=(1-sin((x*%pi)/180))/(1+sin((x*%pi)/180)) //(coefficient of active earth pressure)
 Ka  = 

   0.6557502


--> Pa1=(Ka*Y*z1)-(2*c*sqrt(Ka)) //(stress at top)
 Pa1  = 

  -32.391361


--> Pa2=(Ka*Y*z2)-(2*c*sqrt(Ka)) //(stress at bottom)
 Pa2  = 

   14.822652


--> Zc=(2*c)/(Y*sqrt(Ka)) // (depth of crack)
 Zc  = 

   2.7442159


--> Hc=(4*c)/(Y*sqrt(Ka)) // (maximum depth of unsupported excavation)
 Hc  = 

   5.4884318


--> //results

--> printf('stress at top layer=%f',Pa1) //KN/m^2
stress at top layer=-32.391361
--> printf('stress at bottom layer=%f',Pa2) //KN/m^2
stress at bottom layer=14.822652
--> printf('maximum depth of potential crack=%f',Zc) //m
maximum depth of potential crack=2.744216
--> printf('maximum depth of unsupported excavation=%f',Hc) //m
maximum depth of unsupported excavation=5.488432