The most probable explanation is that the chemical pressure is causing a breaking of symmetry.  The reason, in general, for small lattice changes when you replace La with Eu is that their "sizes" (i.e. atomic radii) are slightly different. When you replace just one atom, you haven't actually reduced the effective symmetry of the system (i.e. the x and y directions are still similar so for example if it were present before the C4z rotation (fourfold rotation about the z axis)  is still present.  However, once you've introduced additional dopant atoms to the system, you are causing a break in the symmetry between these axis.  This causes there to be added chemical pressure in one of the directions.  Beyond what you've already seen, you can also expect the result to be different (lattice constant and total energy) if you place both dopant atoms along the z-axis rather than along the x axis.  Even if you were to randomly place the two dopant atoms in the cell, you will find that some arrangements may or may not break the tetragonality of the cell when the lattice vectors are relaxed.

That being said, it also depends what you are trying to model as to whether or not you should let the lattice constants relax or only the internal coordinates.  For example the cells above represent 6.25%, 12.5%, and 18.75% doping. If you are truly trying to model doping concentrations this high then it is fine to let the lattice constants relax (and at really high percentages you definitely should). However, if you are really interested is lower percentage dopants (i.e. isolated dopants) but the simulation cell size is limited by computational resources then you should keep the lattice constant fixed.

Thank you Dr. Andrew O'Hara for very useful discussion forever  !

Further I want to ask one thing I have to report lattice constants but in some doping concentration it is breaking tetragonal symmetry i.e. a and b are not equal so how will I report that lattice constants.

Other thing how will I know what dopand position will not break the symmetry.

I would just report both a and b separately in a table if tetragonal symmetry were broken. However, it really depends how significantly broken they are. For example, in your second CONTCAR, the difference in a and b is only 0.006 Angstroms; depending upon your convergence criterion this may or may not be a significant difference in a and b.

As for no breaking the tetragonal symmetry, this will always happen when you place certain numbers of dopant atoms in the supercell you have. One thing you could try is for a given percentage of dopant atoms, try several configurations and pick the lowest energy configuration to report. Or if they are similar in energy then some sort of average.