# Solving equation with multi variable matrix input

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]
]
]

• You don't need to discretize the formula yourself, check the document of DensityPlot carefully. Also, notice e.g. E and N are built-in symbol in Mathematica (See their color? They're black, rather than blue ), you can't use them as variable names. Oct 30 '18 at 12:37
• What do you want to determine from these equations? Oct 30 '18 at 13:15
• So I just want to graph the equations. Im not looking for any numerical results. The formulas will tell me the stress and shear below the surface of a material, and I'd like to see the counter plots.
– Alc
Oct 30 '18 at 13:19

After correcting all errors and normalization to A, we have

R = .1;
E1 = 200*10^9;
E2 = 550*10^9;
P = 1000;
V1 = 0.3;
V2 = 0.3;
E3 = 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];

m = Sqrt[.5*(((1 - i^2 + j^2)^2 + 4*i^2*j^2)^(.5) + (1 - i^2 + j^2))];
n = Sqrt[.5*(((1 - i^2 + j^2)^2 + 4*i^2*j^2)^(.5) - (1 - i^2 +
j^2))];
SX = (-P)*m*((1 - ((j^2 + n^2)/(m^2 + n^2))) - 2*j);
SZ = (-P)*m*((1 - ((j^2 + n^2)/(m^2 + n^2)))); SY = V1*(SX + SZ);
TXZ = (-P)*n*((m^2 - j^2)/(m^2 + n^2)); {DensityPlot[
SX/P, {i, -2, 2}, {j, -2, 0},
PlotLabel -> "\!$$\*SubscriptBox[\(\[Sigma]$$, $$x$$]\)/P",
PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}],
DensityPlot[SY/P, {i, -2, 2}, {j, -2, 0},
PlotLabel -> "\!$$\*SubscriptBox[\(\[Sigma]$$, $$y$$]\)/P",
PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}],
DensityPlot[SZ/P, {i, -2, 2}, {j, -2, 0},
PlotLabel -> "\!$$\*SubscriptBox[\(\[Sigma]$$, $$z$$]\)/P",
PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}],
DensityPlot[TXZ/P, {i, -2, 2}, {j, -2, 0},
PlotLabel -> "\!$$\*SubscriptBox[\(\[Tau]$$, $$xz$$]\)/P",
PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}]}
{ContourPlot[SX/P, {i, -2, 2}, {j, -2, 0},
PlotLabel -> "\!$$\*SubscriptBox[\(\[Sigma]$$, $$x$$]\)/P",
PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}],
ContourPlot[SY/P, {i, -2, 2}, {j, -2, 0},
PlotLabel -> "\!$$\*SubscriptBox[\(\[Sigma]$$, $$y$$]\)/P",
PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}],
ContourPlot[SZ/P, {i, -2, 2}, {j, -2, 0},
PlotLabel -> "\!$$\*SubscriptBox[\(\[Sigma]$$, $$z$$]\)/P",
PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}],
ContourPlot[TXZ/P, {i, -2, 2}, {j, -2, 0},
PlotLabel -> "\!$$\*SubscriptBox[\(\[Tau]$$, $$xz$$]\)/P",
PlotLegends -> Automatic, FrameLabel -> {"x/a", "z/a"}]} 