Remove["Global`*"];​​​​Needs["NDSolve`FEM`"]​​​​bmesh=ToBoundaryMesh["Coordinates"{{0.5*10^-3,0.},{1.5*10^-3,0},{20*10^-3.,0},{20*10^-3,7*10^-3},{1.5*10^-3,7*10^-3},{0.5*10^-3,7*10^-3}},"BoundaryElements"{LineElement[{{1,2},{2,5},{5,6},{6,1},{2,3},{3,4},{4,5}}]}]​​bmesh["Wireframe"]​​​​mesh=ToElementMesh[bmesh,"RegionMarker"{{{0.6*10^-3,0.001},1,0.00000015},{{5*10^-3,0.001},2,0.000002}},MaxCellMeasure0.00000051];​​Show[mesh["Wireframe"],mesh["Wireframe"["MeshElementMarkerStyle"Blue]]]​​​​(*PDElösen*)​​lam1=8;​​lam2=20;​​q=695734.466;​​​​boundaries={z(7*10^-3),z0,r0.5*10^-3,r1.5*10^-3,r20*10^-3};​​​​ufun=NDSolveValue[{-D[u[r,z],{r,2}]-D[u[r,z],{z,2}]-1/r*D[u[r,z],{r,1}]NeumannValue[q/lam1,boundaries[[3]]],{DirichletCondition[u[r,z]25,boundaries[[5]]]}},u,{r,z}∈mesh]​​​​Plot3D[ufun[r/10^3,z/10^3],{r,0.5*10^-3*10^3,20*10^-3*10^3},{z,0,7*10^-3*10^3},ColorFunction"TemperatureMap",AxesLabel{"r / mm","z / mm","T / °C"},PlotLabel"Stationäres Temperaturprofil im Kristall (90*0.17 W Heizleistung)"]
Out[]=
ElementMesh[{{0.0005,0.02},{0.,0.007}},Automatic]
Out[]=
Out[]=
Out[]=
InterpolatingFunction
Domain: {{0.0005,0.02},{0.,0.007}}
Output: scalar

Out[]=