The problem is basically that nestedseq was calculating the expression first with a symbolic variable and then evaluating at the values in.felems.  That doesn't work with Determinant[GF4].   Try this one:

> nestedseq:= proc(numtimes::nonnegint, middle::uneval,felems)
  local k,R;
  if numtimes = 0 then R:=eval(middle)
  else
     R:= NULL; 
     for k in felems do
     assign(cat('j',numtimes),k);
     R:= R,nestedseq(numtimes-1, middle,felems);
     end do
  end if;
  unassign(cat('j', numtimes));
  R;
  end proc;

> nestedseq(4, Determinant[GF4](Matrix([[j1,j2],[j3,j4]])),
      [0,1,z,z+1]);

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, z, z, z, z, z+1, z+1, z+1, z+1, 0, 0, 0, 0, z, z, z, z, z+1, z+1, z+1, z+1, 1, 1, 1, 1, 0, 0, 0, 0, z+1, z+1, z+1, z+1, 1, 1, 1, 1, z, z, z, z, 0, 1, z, z+1, 0, 1, z, z+1, 0, 1, z, z+1, 0, 1, z, z+1, 0, 1, z, z+1, 1, 0, z+1, z, z, z+1, 0, 1, z+1, z, 1, 0, 0, 1, z, z+1, z, z+1, 0, 1, z+1, z, 1, 0, 1, 0, z+1, z, 0, 1, z, z+1, z+1, z, 1, 0, 1, 0, z+1, z, z, z+1, 0, 1, 0, z, z+1, 1, 0, z, z+1, 1, 0, z, z+1, 1, 0, z, z+1, 1, 0, z, z+1, 1, 1, z+1, z, 0, z, 0, 1, z+1, z+1, 1, 0, z, 0, z, z+1, 1, z, 0, 1, z+1, z+1, 1, 0, z, 1, z+1, z, 0, 0, z, z+1, 1, z+1, 1, 0, z, 1, z+1, z, 0, z, 0, 1, z+1, 0, z+1, 1, z, 0, z+1, 1, z, 0, z+1, 1, z, 0, z+1, 1, z, 0, z+1, 1, z, 1, z, 0, z+1, z, 1, z+1, 0, z+1, 0, z, 1, 0, z+1, 1, z, z, 1, z+1, 0, z+1, 0, z, 1, 1, z, 0, z+1, 0, z+1, 1, z, z+1, 0, z, 1, 1, z, 0, z+1, z, 1, z+1, 0

Something like this perhaps?

> nestedseq:= proc(numtimes::nonnegint, middle)
    if numtimes = 0 then middle
    else seq(eval(nestedseq(numtimes-1, middle)),
              cat('j',numtimes)=1..3)
    end if
  end proc;

> nestedseq(3, b(j1,j2,j3));

b(1,1,1), b(2,1,1), b(3,1,1), b(1,2,1), b(2,2,1), b(3,2,1), b(1,3,1), b(2,3,1), b(3,3,1), b(1,1,2), b(2,1,2), b(3,1,2), b(1,2,2), b(2,2,2), b(3,2,2), b(1,3,2), b(2,3,2), b(3,3,2), b(1,1,3), b(2,1,3), b(3,1,3), b(1,2,3), b(2,2,3), b(3,2,3), b(1,3,3), b(2,3,3), b(3,3,3)

Normally,  GRADIENT(SK) should not contain any D's unless SK involves functions that haven't been assigned or that Maple doesn't know how to differentiate.  Can you give a specific example (preferably as simple as possible) where you get those D's?

For example, if u(x,t) is one of the dependent variables,

> solu:-plot(u(x,t), t = 0.1);

> eqs:= [(-42.62500000*omega^2+1400.)*A[1]+
     (-12.00000000*omega^2-1000.)*A[2] = 0,
      (-12.00000000*omega^2-1000.)*A[1]+
     (-5.475000000*omega^2+1000.)*A[2] = 0];

> LinearAlgebra[GenerateMatrix](eqs, [A[1],A[2]])[1];

> with(plots): with(plottools):
  display(transform((x,y)->[x,y,0])(disk([0,0],2, colour=red)));

The point of the exercise, I suspect, is to notice just how bad a method of computing Pi this is (it's called Gregory's series).  The error in the N'th partial sum of an alternating series such as this is bounded by the size of the N'th term (and typically not much smaller than that).  To get Pi accurate to within 5*10^(-6), you'd want your series for Pi/4 accurate to within 1.25*10^(-6), so try N = 400000.

> evalhf(4*add((-1)^k/(2*k+1), k = 0 .. 400000));

I assume you mean all those "lembda" to be "lambda".  But the main problem is that your differential equations are undefined when one of the variables is 0, and that's what your boundary conditions are trying to enforce.  In order to use numerical methods, I think you'd have to reformulate your system in some way that removes those singularities.

1) Post text rather than images.

2) Do not use brackets [] for grouping in Maple.  Those are reserved for lists and indices.  Use ordinary parentheses ().

3) Don't forget the multiplication sign *.

4) The exponential function is exp.  e is just another variable.

Your differential equation is a quadratic in y'(x), so for any x and y(x) there will generally be two possible values for y'(x).  Which do you want?  You could do something like this:

> with(DEtools):
  eq:= diff(y(x),x)^2 - 4*x*diff(y(x),x) + 4*y(x) = 0;
  S:= solve(eq, diff(y(x),x));
  P1:= DEplot(diff(y(x),x) = S[1], y(x), x = -2 .. 2, y = -2 .. 2, 
     colour=red):
  P2:= DEplot(diff(y(x),x) = S[2], y(x), x = -2 .. 2, y = -2 .. 2, 
     colour=blue):
  plots[display](P1,P2)

My guess is that this has something to do with the T(t) or m(t).   In a simpler example:

> dsolve({Diff(y(t),t,t)=t,y(5)=0,D(y)(0)=0},y(t),numeric,
     output=array([seq(k,k=0..5)]));

        [                 [         d      ]                 ]
        [                 [t, y(t), -- y(t)]                 ]
        [                 [         dt     ]                 ]
        [                                                    ]
        [[0.    -20.8333333333333252             0.         ]]
        [[                                                  ]]
        [[1.    -20.6666666666666572    0.499999999999999888]]
        [[                                                  ]]
        [[2.    -19.5000000000000000    2.00000000000000044 ]]
        [[                                                  ]]
        [[3.    -16.3333333333333322    4.49999999999999912 ]]
        [[                                                  ]]
        [[4.    -10.1666666666666696    8.00000000000000178 ]]
        [[                                                  ]]
        [[5.             0.             12.5000000000000036 ]]

On the other hand, the same with method=gear produces an error:

Error, (in dsolve/numeric/type_check) boundary value problem 
cannot be solved using the gear method

The point being that the gear method is only for initial value problems, not for boundary value problems. 

The way to specify an initial condition for a derivative is using D rather than diff, e.g.

D(y)(0) = 0

There's no stochastic differential equation here, just a partial differential equation.  Using Fourier transforms is one of the standard ways of solving this.
 

Something like this, perhaps:

> # define the function f
   f:= (x,y) -> x/10 + y^2;   
   # compute your points
   for i from 1 to 5 do
     X := i * 0.1 - 0.3;   
     for j from 1 to 10 do
       Y:= j * 0.1 - 0.5;
       Z:= f(X,Y);
       Data[i,j]:= [X, Y, Z];
     end do
   end do:
   # plot the data
   with(plots):  
   surfdata([seq([seq(Data[i,j],j=1..5)],i=1..5)],
      axes=box, labels = ["x","y","z"]);

A somewhat more efficient way, when the x and y values are on a regular grid:

> # define the function x
  f:= (x,y) -> x/10 + y^2;
  # compute the z values 
  ZData := Array(1..5, 1..10, 
             (i,j) -> f(i*0.1 - 0.3, j*0.1 - 0.5),
            datatype=float[8]):
  # plot the data
  with(plots):
  surfdata(ZData, -0.2 .. 0.2, -0.4 .. 0.5, 
     axes = box, labels = ["x","y","z"]);

  Or in such a simple example, you could just use plot3d:

> plot3d(f(x,y), x = -0.2 .. 0.2, y = -0.4 .. 0.5, 
     grid = [5, 10], axes = box, labels = ["x","y","z"]);

Or, if f is a function that doesn't work for symbolic arguments,

> plot3d(f, -0.2 .. 0.2, -0.4 .. 0.5, 
     grid = [5,10], axes=box, labels = ["x","y","z"]);

Do you mean something like this?

> with(plots):
  Data:= [[1,2,3,5],[4,5,6,6],[7,8,8,7],[9,8,7,6],[7,5,3,1]]:
  surfdata(Data, axes=boxed, labels=[x,y,z]);