I am trying to solve this system of elasticity partial differential equations in 2D in cylindrical coordinates, where the model is axysimmetric, around the $z$ axis.
$$\frac{1}{r} \frac{\partial (r\sigma_{rr})}{\partial r} +\frac{\partial \sigma_{rz}}{\partial z} =f_{r}$$
$$\frac{1}{r} \frac{\partial (r\sigma_{rz})}{\partial r}+\frac{\partial\sigma_{zz}}{\partial z}=f_{z}$$
Let us consider a solution of eq. for $f_{r}=0$ and $f_{z}=-\rho g$ with the property that at the top plane $z=h$: $\sigma_{zz}=-\frac{F}{A}$. If $A(z)$ denotes the area of horizontal cross-section in $z$-direction. $F=F(z)$ is the magnitude of force acting uniformly at the cross section in $z$-direction.
System should be solved for the stresses over entire domain: $\sigma_{rz}(r,z), \sigma_{zz}(r,z), \sigma_{rr}(r,z)$. The boundary conditions are:
The cylindrical surface at $r=a$ and the upper plane at $z=h$ are traction free:
$\left(\begin{matrix}\sigma_{rr}\big|_{r=a}\\ \sigma_{rz}\big|_{z=h}\\ \sigma_{rz}\big|_{r=a}\\\sigma_{zz}\big|_{z=h}\end{matrix}\right) =\left(\begin{matrix}0\\ 0\\ 0\\-\frac{F_0}{A_0}\end{matrix}\right)$
For the bottom plane smoothly contact with the rigid base $z=0$.
$\left(\begin{matrix}\sigma_{rz}\big|_{z=0}\\ \sigma_{zz}\big|_{z=0}\\ \end{matrix}\right)
=\left(\begin{matrix}0\\ 0\\ \end{matrix}\right)$
I tried with Maple, but seems that Maple does not solve a system of elliptic PDEs, only single elliptic PDE. 
Then I tried with Mathematica, by using the following code:
k = 7800*9.8;
pde1 = D[r*Subscript[\[Sigma], rr][z, r], r]/r +
D[Subscript[\[Sigma], rz][z, r], z] == 0;
pde2 = D[r*Subscript[\[Sigma], rz][z, r], r]/r +
D[Subscript[\[Sigma], zz][z, r], z] == -k;
constrains = {Subscript[\[Sigma], rr][z, 10] == 0,
Subscript[\[Sigma], rz][10, r] == 0,
Subscript[\[Sigma], rz][z, 10] == 0,
Subscript[\[Sigma], rz][0, r] == 0,
Subscript[\[Sigma], zz][10, r] == 1000, Subscript[\[Sigma], zz][0, r] == 0};
domain = ImplicitRegion[0 <= z <= 10 && 0 <= r <= 10, {z, r}];
sol = NDSolve[{pde1, pde2, constrains}, {Subscript[\[Sigma], rr][z,r], Subscript[\[Sigma], zz][z,r], Subscript[\[Sigma], rz][z, r]}, {z, r} \[Element] domain]
And i get the answer: NDSolve::underdet: There are more dependent variables, {Subscript[\[Sigma], rr][z,r],Subscript[\[Sigma], rz][z,r],Subscript[\[Sigma], zz][z,r]}, than equations, so the system is underdetermined.
Can this be solved, because I have a system with 3 unknowns and 2 equation, my system is undetermined, or if i can add some more condition?
Isn't the solution the following?$$\sigma_{zz}=-\frac{F}{A}+(h-z)\rho g$$ $$\sigma_{rr}=\sigma_{rz}=0$$