You want to use Matrices instead of subscripted names. It'll be much faster. Then your code will look like this:

restart:
dellT:= 0:  h:= 1:
tile:= Matrix(150$2, datatype= float[8]):
#You'll need to fill the above matrix with initial values.
tile2:= Matrix(tile):
to 20 do    
     for j from 2 to 149 do
          for k from 2 to 149 do  
               tile2[j,k]:= dellT/h^2*(tile[j+1,k]-4*tile[j,k]+tile[j-1,k]+tile[j,k-1])+tile[j,k]  
          end do
     end do;  
     tile:= copy(tile2)
end do:

Assuming that you'll be generating a lot of these pairs, you'll want to avoid regenerating the generating procedure, which is inefficient. This can be done with a module.

GenXY:= module()
local
     R:= RandomTools:-Generate(list(integer(range= 1..7), 2), makeproc),
     ModuleApply:= proc()
     local x:= 0, y:= 0;
          while x=y do  (x,y):= R()[]  end do;
          (min(x,y), max(x,y))
     end proc
;
end module:
     
(x,y):= GenXY();

                              2, 7

This Answer is essentially the same as John's, but uses a Vector of lists instead of a list of lists. You could also use a Vector of Vectors or a Matrix. The best format to store the data should be determined by how many times you do these shifting operations. Shifting operations on lists are inefficient.

block:= < [0,0,1,1,0,0,1],[0,0,1,1,1,0,0],[0,1,0,1,0,1,0],[1,0,0,1,1,1,0]>;

block:= ArrayTools:-CircularShift(block, 1);

plots:-display(
     plots:-pointplot3d([[0,0,12]], color= red, symbolsize= 24),
     plots:-spacecurve([4*t,-2*t,2*t], t= 0..10, thickness= 4),
     axes= boxed
);

So, is f3 the result (or one of the results) of a dsolve(..., numeric) computation? And you want to compute the integrals a31 and a32 after doing the dsolve, right? Shouldn't that f3(x) be f3(theta)?

Assuming that the answers to those questions are yes and that you've correctly extracted the f3 from the dsolve solution, then you compute the second integral with

a32:= evalf(Int(unapply(chi*g3/(1-f3(theta)*g3)^4, theta), a..1));

And likewise for the first integral.

@mskalsi

[Edit: This is an Answer to the OP's followup question posted as a Reply to my first Answer below.]

Yes, that can be done in a loop, like this:

Download Loopwise.mw

You need to add the word identity:

solve(identity(f(t), t), {a,b});

Better than a loop, you can get a tabular layout by putting the values in a Matrix, which is trivial:

Matrix(5, 5, F);

By the way, I don't see anything recursive about your Question.

I assume in the code below that x and y are unassigned names and that X, Y, and Z are Vectors containing your data.

V:= [x,y]:  deg:= [3,3]:
poly33:= [op(expand(mul(add(V[j]^k, k= 0..deg[j]), j= 1..nops(V))))]:
fitresult:= Statistics:-LinearFit(poly33, <X|Y>, Z, V);

Use *~:

with(Units[Standard]):
X:= [1,6,2]*~Unit(m):
Y:= [2,1,4]*~Unit(Hz):
X*~Y;

Let ex be the expression. Then do

C:= coeffs(expand(ex), indets(ex, indexed), 'T'):
eval(<U[1],U[2]>.<a[1] | a[2]>, [T=~C]);

Although you can often get away with using some other of Maple's numerous derivative notations, I recommend that initial and boundary conditions be specified in D notation. The following command returns the particular solution without complaint:

dsolve({eq, y(0)=a, D(y)(0)=0, (D@@2)(y)(0)=0}, y(x));

I put your code into a module in such a way that it handles all the issues that you raised.

koch:= module()
local
     points, npts, t,
     Init:= proc(n::nonnegint)
          (points, npts, t):= (table([1= [0,0]]), 1, [1/3^n,0]);
          [][]
     end proc,

     forward:= proc()
          local elem:= points[npts];
          npts:= npts+1;
          points[npts]:= elem +~ t;
          [][]
     end proc,
   
     turnleft:= proc() t:= [-t[2], t[1]]; [][] end proc,
     turnright:= proc() t:= [t[2], -t[1]]; [][] end proc,
     turndegree:= proc(deg::numeric)
     local c,s;
          (c,s):= evalf([cos,sin](deg*Pi/180))[];
          t:= [t[1]*c-t[2]*s, t[1]*s+t[2]*c];
          [][]
     end proc,

     Rand:= rand(2),

     Recurse:= proc(n)
          local r;
          if n=0 then return forward() end if;
          r:= `if`(Rand()=0, -1, 1);
          thisproc(n-1), turndegree(r*60),
               thisproc(n-1), turndegree(-r*120),
               thisproc(n-1), turndegree(r*60)
     end proc,   

     ModuleApply:= proc(n::nonnegint)
          Init(n);
          Recurse(n);
          plot(
               convert(points,list),
               thickness= 2, scaling= constrained, axes= none, _rest
          )
     end proc
;
end module:

Now to use it, you simply need to enter, for example,

koch(5);

You can also include additional plot options, for example,

koch(4, color= violet, thickness= 3);

length(tekssifer);

All operations can be done with a one-line procedure if you're willing to input the operations in prefix form. The previous two Answers also require prefix form anyway, but my prefix form is identical to the original operations.

`&B`:= (ex::uneval)-> convert(eval(map(convert, ex, decimal, 2)), binary):

Your examples:

&B `+`(1101, 111);

     10100

&B `-`(11000, 1011);

     1101

&B `*`(1011, 1101);

     10001111

Integer division is done in Maple with the command iquo.

&B iquo(10010011, 1011);

     1101