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

Here is your program with corrected syntax and plot of U:

restart:
H:= 0.5:
q__0:= 10^5:
Es:= 4*10^9:
p__0:= 10^6:
q__s:= 10^10:
C:= q-> q^(H+1.5):
G:= q-> int(qs^3*C(q), qs= q__0..q):
P:= q-> 1/sqrt(G(q)):
p:= xi-> p__0/P(xi*q__0):
w:= (q,xi)-> 1/(Pi*(int(C(qs)*qs^3, qs= xi*q__0..q)))^(1/2):

U:= xi-> int(q*C(q)*w(q,xi)*Int(exp(-(w(q,xi)*ps/Es)^2)/ps, ps= p(xi)..infinity), q= xi*q__0..q__s):

plot(U, 1.5..5);

How do you expect the plotting routines to deal with the imaginary parts? If you use Re to extract the real part of your function, then you can get a plot. Like this:

ex1:= (x,t,z)-> Re(-1.132e11*exp(9.9e6*x + I*(1.95e6*z-2.98e15*t)));
plots:-implicitplot3d(
     ex1(x,t,z), x= -10..0, t= 0..10, z= 0..10,
     axes= boxed, style= patchcontour, scaling= constrained, shading= z,
     grid= [20$3]
);

PS: Because your function can never be 0, as Preben pointed out, the plot produced by the code above is just random noise.