This is your worksheet with some minor corrections, comments, and three different ways of arriving at the same result (note you didn't say what is the result you expect)

I see you are using primed variables, which is perfectly fine, just recalling the prime in this context does not mean derivative.

Here is the first adjustment: you use the assignment operator :=, but on the left you put functionality: in doing so, you create a function, as you see in the output of your post, but after that, you only use the primed variables as symbols, not as functions, so your definition is not used (this is not related to Physics, but to how computer algebra systems work).

I am then removing the functionality on the lhs so that the primed variables have the value you indicated

You see how the value is now taken into account:

Next you assigned to a sum, I imagine to want automatic simplification ... it is not how it works: if you want the following simplification you need to invoke it. I am changing that := by  = and show

Next you use , not wrong but perhaps more clear: you could have used its value

Here again you see the value of the primed variables, not themselves

Again, the output involves the value of the primed variables

Now the Sum over all the repeated indices

Is this a scalar? Yes, in the sense that there are no tensors around. Is this the result you were expecting? I don't know, it would help if you type in the result you were expecting.

Also, you see that (17) above depends on  and . You can determine the conditions on  such that  depends only on r by differentiating

Any  that cancels both equations result in  depending only on r.

Alternatively, you can compute the same but keeping the primed variables around till the last step, as follows

So do not assign the primed variables. Instead, create a set of substitution equations (i.e. use =, not :=) to be used only when you want

Proceed now using the primed variables themselves instead of their expression in terms of spherical coordinates

You see now the definition in terms of the primed variables

which all depend on the spherical coordinates

Next, the following is an identity, I will assign it here to some name, eq, in order to use it for simplification purposes later

Define now A

Here again you see the primed variables themselves, not their value

All primed ones are functions of the spherical coordinates:

Define now your F

Again, the output involves the primed variables, and their derivatives, denoted with an index

Now the Sum over all the repeated indices

Simplify the above taking into account eq

Express now the primed variables as functions of the spherical coordinates

So your  now is given by

This is the same result obtained following the previous approach, where

A last verification, also resulting in the same: compute with everything inert, to revise by eye each step. This is the Lagrangian

You see that the su2 indices, last index in F, run from 1 to 3 only

Compute the value of each of these inert F

Now for the contravariant ones

At this point, you have the Lagrangian in terms of inert F in (47), so that you can follow by eye that the computation is running as expected, and in (49),(50) you have the value of each inert F. So go ahead and de-inertize all F in (47)

Introduce the value of the primed variables in terms of spherical coordinates to get your

This is again equation (17), the Lagrangian you got in the first approach

You know, PDE solutions depend on arbitrary functions, so specializing them you have infinitely many particular solutions. The one you show is very so, it has no arbitrary functions and no arbitrary constants. The pdsolve command attempts to find general solutions, and if that is not possible, then complete solutions or else it does not attempt to find all and any possible particular solution, but for separable (simple or more general forms) or travelling wave ones.

Anyway, for more specific particular forms of solutions, e.g. polynomial ones as the one you show, there are specific commands in PDEtools, in this case PDEtools:-PolynomialSolutions.

So here are some more general - still particular - solutions, these ones depending on arbitrary constants.

This one depends on 2 arbitrary constants

This other one depends on 4 arbitrary constants,

This other one depends on 6 arbitrary constants; I omit the output to avoid cluttering.

The 6 arbitrary constants (you could also use suffixed(c__))

There are different things you may want regarding the covariant derivative of Einstein's tensor. Below you have a sample of them.

Note that after a call to CompactDisplay , derivatives are displayed indexed; to see what is behind this display, use show

The covariant derivative, you may want to compute with it symbolically (simplification of expressions involving it is invoked using Physics:-Simplify).

(7) has three indices, the default display is a slice

This is practical

This gives you a slice display

This allows you to explore each slice, having  or  varying from 1 to 4

For details, or in general, check the help pages for Physics:-D_ and Physics,Tensors.

Note that r is a coordinate, do not define it also as a tensor. Use another letter, for instance R

Also, define all the components in one go instead of one component at a time

Check it out

The same for , define all its components in one go, it is simpler, so comment these lines, then define it

Here is the above in one go: use LeviCivita with four indices, as in

Check it out

Your next definition: my suggestion to you is again to use  with four indices, the forth being the number 4

Next you could use local D but you will find that the results of your computation involve the D Maple differentiation operator (check its help page), so to avoid this problem use _D, and make it point to any convenient letter you can get from the palettes (the column at our left), for instance

Now use _D all around for your covariant derivative operator

To see the components of any tensor, you can use , or directly this

At this point, your question is about computing the functional derivative with respect to . The help page for Fundiff says that you need to specify a different variable for the differentiation function, and then you get  the expected Dirac function, as in

So not the same variable for the function - you would receive this kind of telling error message

Also, you do not need to put the Lagrangian density within an integral in order to compute a functional derivative: Fundiff and Int commute, so you can compute, directly

Rewrite the D derivatives using diff

Now you can now integrate this result with any limits of integration that you prefer. Note however that you have two Dirac functions

and no Dirac function for the other variables, so how do you intend to integrate in the other variables? Also, are these Dirac functions what you expect?  The computer can perform for you a computation that you know how to do it. That is important also to know whether the result is actually what you wanted to compute.

I suggest you to give a look for instance at at this link with brief information on "The Dirac Delta Function in Three Dimensions", or search the web for related material. Using Intc as you did will integrate over these coordinates as if they were Cartesian, from - to + infinity. If after formulating the problem in a way that you are sure you are computing what you want to compute you still have a doubt, feel free to post here again and we take if from there.
 

(Please remember to post a worksheet instead of text and screen shots.)

There is no bug. This is the answer to your question, and at the end two suggestions to avoid these subtleties.

From your signature, you are choosing position 1 to be timelike, and 2, 3 and 4 to represent the three spacelike components. With that in mind, take a look at the components of

So the Pauli matrices are ,  and , as you can also see in this output

So what you do expect for the components of  where Psigma is indexed with a space index a?

The above turns on the lights regarding the shifted value of the space index a : the tensor you are indexing,  is a 4-dimensional tensor, and being that the 1st component is timelike (not Physics's default but your choice of signature), the space components are given by shifting the 4-dimensional index by one.

How about a true space tensor? Let's define one of them, another mixed, and one all spacetime for experimentation

For , a space tensor, the components are given by

You see they are not shifted. The index a runs from 1 to 3, not from 1 to 4 such that you are only interested in the components from 2 to 4.

How about the mixed tensor ?

Again OK, and you see: there are 4 lines indexed from 1 to 4, and 3 colums indexed from 1 to 3, because a is a space index. How about , where in  we want only the space part of the spactime index ?

This is again correct, but the 3 columns are indexed from 2 to 4, not from 1 to 3 as in (8) for . I hope it is clear now: when you set a signature with timelike component in position 1, the space part of a spacetime index is given by shifting the value by 1 to skip the timelike component in position 1.

With these things in mind lets give a closer look at your example where you define a covariant derivative operator as . The first problem is that you use the letter D, which is a Maple command representing differentiation. If you proceed as you did, you effectively not have the D command anymore, which is bad. Instead, I suggest you first input

Now you can use D without cancelling the Maple differentiation command D

Check

Let's sum over the repeated index a. In view of the comments lines above, what do we expect? In , the index a is a space index, running from 1 to 3. In , the index a represents the space part of a spacetime index so it runs from 2 to 4. Then

That is the result you are getting when entering and that motivated your question: the index a in  runs from 1 to 3 but the index a in  runs from 2 to 4.

Of course you could have all this in mind, and the computations will run correctly, you don't need to worry. Or the details are not in your mind and get in doubt (your post). I have some suggestions to avoid falling in these subtleties

Compare. These are the current components for your A defined as a mixed components tensor

with these other ones after redefining A with two spacetime indices

I believe the above is what you wanted. So now you have

and you get what you where expecting  for .