Thank you for your help. 
It is not exactly what I would expect. I would expect my Lagrangian density to depend only on r, due to the radial symmetry of my problem. But that's up to me, I made a small mistake in the definition. Unfortunately I only know the final result of the problem (plot of the equation of motion, after solving the Euler lagrange equations) but with your help I managed to get a system of equations, which is very good! 

One more question about the difference between := and =.  If I use r := sqrt(x^2+y^2+z^2) then d/dx r = x/sqrt(x^2+y^2+z^2) but if I say r = sqrt(x^2+y^2+z^2) then d/dx r = 0, what exactly is the difference?

Many thanks for the detailed help! 

I wouldn't expect a Dirac function, that's right, I probably made a mistake somewhere. The Lagrangian density is also not yet complete, my Higgs field is still missing. I just tried to do a variational calculation to see how it works in general.  

Thank you very much! 

Thank you, you have helped me a lot! 

@ecterrab 

Thank you for the quick reply. I have only been using maple for a very short time. Do you mean this page with the physics help page? 

https://www.maplesoft.com/support/help/Maple/view.aspx?path=Physics/Tensors&cid=412