Looks like you have some hybrid Mathematica/Maple syntax. Maple uses round parentheses instead of square brackets. The command for extracting the real part of an expression is a combination of evalc and Re:

evalc(Re(sqrt(x+I*y)));

The evalc tells Maple to assume that all variables are real. The Re alone wouldn't return the above answer because that assumption wouldn't have yet been made.

Here is the algorithm coded in Maple and a test instance. Any advice on how I can do more of the Matrix operations inplace would be much appreciated.

Download Fuzzy_clusters.mw

You want

plot(eval(C, k= k/K), k= 0*K..4*K, labels= ['K', 'C']);

You can't have more than one independent variable in an ODE system that is solved by dsolve. You have x and y as independent variables.

Use the bits option of Bits:-Split. For example,

mb:= Bits:-Split~(convert(message, bytes), bits= 8);

You need a multiplication sign:

g:= x-> exp(x^2*(x-a));

This is done by computing the exclusive-or (xor) of the two bit strings and counting the ones in the result.

CountOnes:= proc(n::nonnegint)
local q:= n, c:= 0;
     while q > 0 do  c:= c+irem(q,2,'q') end do;
     c
end proc:
 
CountFlips:= (X::{list, Matrix, Vector}, Y::{list, Matrix, Vector})->
     CountOnes(Bits:-Xor((Bits:-Join@convert)~([X,Y], list)[]))
:

CountFlips(<1,0,0,0,0,0,0,0>, <1,1,1,1,0,0,1,1>);

     5

Change your definition of B_interp to

B_interp:= (sr, e_g)->
     CurveFitting:-ArrayInterpolation(
          [SR, E_G], Vert_Coef, [[sr], [e_g]], method= linear
     )[1,1]
:

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});