In the usual sense (Lebesgue/improper Riemann or what have you) the function exp(I*(k-q)*x) is not integrable on the interval -infinity..infinity (of course).

Thus it is entirely appropriate for Maple to return 'undefined'.

However, FunctionAdvisor is aware of the connection to Dirac:

FunctionAdvisor(Dirac,definition);
FunctionAdvisor(Dirac,integral_form);

f[1] := x[2]+x[3]+x[5]+x[6]=16; f[2] := x[1]+x[2]-x[3]+x[4]+3*x[6]+x[7]=28;
var:=indets({f[1],f[2]},And(name,Not(constant)));

#The maximal index:
max(op(op~(var)));
#The corresponding variable:
op(select(has, var,%));

You can use implicitplot3d.

with(plots):
f:=x^a * y^b * z^c ;
fk:=eval(f,{a=1/2,b=1/2,c=1/2});
implicitplot3d(fk=50, x=0..40 ,  y=0..40, z=0..40,axes=boxed);
implicitplot3d([fk=50,fk=100,fk=150] ,x=0..40 ,  y=0..40, z=0..40,axes=boxed);
animate(implicitplot3d,[fk=L, x=0..40 ,  y=0..40, z=0..40,axes=boxed],L=50..150);
animate(implicitplot3d,[eval(f,{a=1/2,b=1/2})=50, x=0..40 ,  y=0..40, z=0..40,axes=boxed],c=0..1);

Without redirecting I tried your code and it produced the attached matlab output together with the warnings

Warning, the function names {`theta[1]`, `theta[2]`, `theta[3]`, `theta[4]`} are not recognized in the target language
Warning, the following variable name replacements were made: ["cg0", "cg2", "cg4"] = ["Δx", "Δy", "Δz"]

matlab.txt

Have you looked in the help pages associated with dsolve/numeric starting with

?dsolve,numeric

I'm not sure that the following is what you are thinking about, but you could try it:

plots:-animate(plot,[[x^2+y,x*y,x=-1..1]],y=-1..1);
plots:-animate(plot,[[x^2+y,x*y,x=-1..1]],y=-1..1,trace=24);

sol := solve({x+y=a, x-y=b}, {x,y});

subs(sol,x/y);
#Alternatively
eval(x/y,sol);


Adding an element to a list could be done as shown:

c:= [1,2,4,4,-3];
[op(c),Pi];
[op(1..2,c),Pi,op(3..-1,c)];

A new list is created in each case.

Tables are useful in some cases where the size is not determined beforehand:

c:= [1,2,4,4,-3];
T:=table(c);
T[6]:=Pi;
eval(T);
#Inserting an element meant to be between the elements with indices 3 and 4:
T[3.5]:=hh;
sort(op(eval(T)),(x,y)->lhs(x)<lhs(y));
rhs~(%);

Here are two ways:

X:=[seq(i, i=1..8)];

t:=[2,5,6];

#Method 1:

Y:=X:
for i from 1 to nops(t) do Y:=subs(t[i]=NULL,Y) end do;

#Method 2

eval(subs(`=`~(t,nULL),X),nULL=NULL);

If in the line A:=B you intend to make a new matrix called A and having the entries which B has at that time, then you should use A:=Copy(B) instead.

Avoid assigning to theta[eta] at the end, this will implicitly create a table called theta. Call the result something else. The problem is that you are also using theta as the name of the unknown function.

There will be other problems but that one is critical.

Surely you will have to solve numerically, and I suggest replacing the number 26 by a more modest value to begin with (e.g. 6). Also replace infinity in the boundary conditions with a finite but large value like 1000 (or even less).

Something like this will then work

res:=dsolve({ics, ODE1},numeric,abserr=1e-1,maxmesh=8192);

plots:-odeplot(res,0..1000);

You can use coeffs.

The following lengthy version does not use that the unknowns are indexed except when defining X, at which time you decide their ordering.

restart;

z:= 7*x[1]+11*x[3]-10*x[4]+26*x[6];

X:=[seq(x[i],i=1..6)];
R:=coeffs(z,X,'t');
t;
T:=table(`=`~([t],[R]));

#As an alternative to the elementwise operation ~ you can use zip:
#T:=table(zip(`=`,[t],[R]));
for u in {op(X)} minus {t} do T[u]:=0 end do;
eval(T);
eval(X,op(eval(T)));

Here is a version which makes u[i], DI1, and DI2 functions, and thereby avoids subs.

Also I have made use of `assuming` and removed evalf.

restart;
lambda:=-1; alpha:=1/2;  N:=2;  m:=ceil(alpha);

                             
u[0]:=t->0;      #initial conditions
for i from 0 to N do
 DI1:=unapply(simplify(1/GAMMA(m-alpha)*int((x-t)^(m-alpha-1)*diff(u[i](t),[t$m]),t=0..x)),x) assuming x>0;
 DI2:=unapply(DI1(t)-t^(1/2)/Pi^(1/2)*(u[i](t)+1)^2,t);
 u[i+1]:=unapply(u[i](x)+simplify(int(1/GAMMA(alpha)*lambda*(x-t)^(alpha-1)*DI2(t),t=0..x)),x) assuming x>0;
end do;

Matrix([seq([ n , n^2, n!, n^2 -n!], n= 1..5)]);
latex(%);

Quoting from the help page for solve,

"If the output of the solve command is a piecewise-defined expression, then the assuming command can be used to isolate the desired solution(s). If the output is not piecewise-defined, in particular, if the output is constant, assumptions on the independent variables may be ignored. If there are parameters in the input equations, the solve command will use those assumptions in its computations."

Try this instead:

solve( {x^2=4,x>0},x);