Hi Barkpuz. I think your code can be implemented in Mathematica. I've coded much, much more complex numerical PDE problems. However, I don't have the time to study your code in detail. I've made some attempt to convert it. I get array over-runs though. I'm doing something wrong. Perhaps you can further improve it. I know, it looks messy. Try and work with it though as this will help you learn Mathematica. Here's some notes:

1. Use Ceiling[arg] for the ceiling function.

2. Arrays use double brackets as in x[[m,n]].

3. Points are constucted as Point[{x,y,z}]

4. To code a colored point list, need the set of points in french brackets with a color preceding the point structure as in:

plist={{Blue,Point[{x1,y1,z1}]},{Red,Point[{x,y,z}]},...}

Then can use Show[Graphics3D[plist]] to show them

5. I converted the code below to "raw input form" in order to cut and paste it. If you wish, you can cut and paste it into Mathematica, select the cell, then choose Cell/Convert To/ Standard Form to recover most of the math format.

6. In Mathematica, never use upper case letters for user variables; may conflict with standard definitions. So I changed A to a.

Code:

s = 20;
tau := 3
sigma := 2
x = Table[1, {m, 1, s + tau + 1}, {n, 1, s + sigma + 1}];
For[i = 1, i <= (1/2)*s*(s + 1), i++,
u = Ceiling[(1/2)*(Sqrt[8*i + 1] - 1)];
m = (1/2)*u*(u + 1) - i; n = u - m;
a = If[m + tau + 1 == n + sigma, 1,
1/((m + tau + 1)*(n + sigma)); ];
x[[m + tau + 1,n + sigma + 1]] =
x[[m + tau + 1,n + sigma]] - x[[m + tau + 2,n + sigma]] -
a*x[[m + 1,n]]; ]
plist = Table[If[m + n <= s + tau + sigma + 2,
If[x[[m,n]] <= 0, {Green, Point[{m - tau - 1,
n - sigma - 1, x[[m,n]]}]},
{Blue, Point[{m - tau - 1, n - sigma - 1, x[[m,n]]}]}]],
{m, 1, s + tau + 1}, {n, 1, s + sigma + 1}]
Show[Graphics3D[plist, PlotRange -> {{-25, 25}, {-25, 25},
{-5, 5}}]]