Background, here are the equations that I am trying to solve:

Where R, E1, E2, V1, V2, P are all user inputs. X/A goes from -2 to 2 and Z/A goes from 0 to -2. Below is the code that I have so far. I created a list of inputs. Then created two arrays for the x and z inputs. The last is where I am having trouble. I’m trying to create a code such that it will hold a value for X constant in SX, SZ, and TXZ and plug in all the values for Z. Then move to the next value for X and plug all the values in for all the Z. The end goal is to create a density plot that for SX, SZ, and TXZ. Thank you!

`R = .1; E1 = 200*10^9; E2 = 550*10^9; P=1000; V1 = 0.3; V2 = 0.3; E = 1/(((1-(V1^2))/E1)+((1-(V2^2))/E2)); A = ((.75*P*R)/(1.61172*10^11))^(1/3); X = Range[-2 A, 2 A, 0.01*3*A]; Z = Range[0,-2 A, 0.005*3*A]; ZZ = ConstantArray[Z[[Range[Length[Z]]]], Length[X]]; XX = ConstantArray[X[[Range[Length[X]]]],Length[Z]]; For[i=1,i=Length[XX], For[j=1,j = Length[ZZ], M = Sqrt(.5*(((A^2-i^2+j^2)^2+4*i^2*j^2)^(.5)+(A^2-i^2+j^2))) N = Sqrt(.5*(((A^2-i^2+j^2)^2+4*i^2*j^2)^(.5)-(A^2-i^2+j^2))) SX = (-P/A)*M*((1-((j^2+N^2)/(M^2+N^2)))-2*N) SZ = (-P/A)*M*((1-((j^2+N^2)/(M^2+N^2)))) SY = V1*(SX+SZ) TXZ = (-P/A)*N*((M^2-j^2)/(M^2+N^2)), DensityPlot[SX/P,XX/A,ZZ/A] ] ] `