Yeah you are wright.

I thought that at least I put an image of an flow chart.

Here is a code written in Mathematica:

**Code:**

state[\[Sigma]10_, \[Sigma]20_, e0_, n0_, c0_, f0_, step_] :=

Module[{matD, F, c, \[Phi], \[Lambda], I1, J2, \[Sigma]1, \[Sigma]2,

Em, \[Nu], matDep, \[Sigma], edf, avek},

c = c0; \[Phi] = f0; \[Sigma]1 = \[Sigma]10; \[Sigma]2 = \[Sigma]20;

Em = e0; \[Nu] = n0;

edf = Table[{If[i == 1, 0, 0.00001], 0}, {i, 1, step}];

\[Sigma] = {{\[Sigma]1}, {\[Sigma]2}};

matD = Em/((1 + \[Nu]) (1 - 2 \[Nu])) ({

{1 - \[Nu], \[Nu]},

{\[Nu], 1 - \[Nu]}

});

Do[\[Sigma] = \[Sigma] + edf[[i]].matD;

F = -((6 c Cos[\[Phi]])/(Sqrt[3] (Sin[\[Phi]] + 3))) + Sqrt[

1/6 ((\[Sigma][[1, 1]] - \[Sigma][[2, 1]])^2 + \[Sigma][[1,

1]]^2 + \[Sigma][[2, 1]]^2)] + ((\[Sigma][[1,

1]] + \[Sigma][[2, 1]]) (2 Sin[\[Phi]]))/(

Sqrt[3] (Sin[\[Phi]] + 3));

term = If[F <= 0, "Elastic", "Plastic"];

avek = If[\[Sigma][[1, 1]] == \[Sigma][[2, 1]] == 0, 0, ({

{(2 \[Sigma][[1, 1]] +

2 (\[Sigma][[1, 1]] - \[Sigma][[2, 1]]))/(

2 Sqrt[6]

Sqrt[\[Sigma][[1,

1]]^2 + (\[Sigma][[1, 1]] - \[Sigma][[2,

1]])^2 + \[Sigma][[2, 1]]^2]) + (2 Sin[\[Phi]])/(

Sqrt[3] (3 + Sin[\[Phi]]))},

{(-2 (\[Sigma][[1, 1]] - \[Sigma][[2, 1]]) +

2 \[Sigma][[2, 1]])/(

2 Sqrt[6]

Sqrt[\[Sigma][[1,

1]]^2 + (\[Sigma][[1, 1]] - \[Sigma][[2,

1]])^2 + \[Sigma][[2, 1]]^2]) + (2 Sin[\[Phi]])/(

Sqrt[3] (3 + Sin[\[Phi]]))}

})];

\[Lambda] = If[F > 0, F/Flatten[Transpose[avek].matD.avek][[1]], 0];

matDep = If[F > 0, matD (( {

{1, 0},

{0, 1}

} ) - avek.Transpose[avek].matD/

Flatten[Transpose[avek].matD.avek][[1]]), matD];

Print[{\[Sigma] // MatrixForm, F, term, \[Lambda],

matDep // MatrixForm}], {i, 1, step}];

]

And here is an image of flow chart:

http://www.imageupload.org/?d=46B3038B1With best regards

freequo