Hello Thomas,
Sorry, I should have been more clear.
You need to pass all the B[i,j] used by your 'W' routine explicitly as
individual parameters.
So if N was 3, you would need to pass:
i,k,j,l,B[1,1],B[1,2],B[1,3],B[2,1],....
or the 4 integer parameters, and the 9 B[i,j] parameters that form the system.
When the numerical DE solver needs to evaluate the equations, it does so
by passing in an array of strictly numerical values, and getting out an
array of strictly numerical values - it cannot handle passing of arrays of
the dependent variables to other procedures.
Now you can explicitly code a general procedure to handle this type of problem
by writing Maple code for the evaluation of the rhs of the DE system.
In your case, here is what I would do (I will generally provide examples with N=3). I have attached the worksheet content below:

1) Determine an ordering of the B[i,j] that maps these to a single sequence

of variables - I would just use the same sequence order that you have

provided with your double sequence call, essentially row-major order.

**> N := 3:**

**> dvars := [seq(seq(B[i,k](t),k=1..N),i=1..N)];**

2) Now code a dsolve/numeric evaluation procedure explicitly to do

what you want:

**> dproc := proc(n,t,Y,YP)**

**> local N,W,B,i,j,k,l;**

**> **

**> # System is n=N^2, so determine N from n**

**> for N from 1 while N^2**

**> if N^2<>n then**

**> error "unexpected number of variables";**

**> end if;**

**> **

**> # Now transfer values to the B[i,j] for convenience**

**> B := hfarray(1..N,1..N):**

**> for i to N do**

**> for j to N do**

**> B[i,j] := Y[N*(i-1)+j];**

**> end do;**

**> end do;**

**> **

**> # Now declare and compute 'W'**

**> W := hfarray(1..N,1..N,1..N,1..N);**

**> # Some complicated computation for W[i,j,k,l] here that can depend on**

**> # the B[i,j], on 't' on N, etc.**

**> # For now, use your example, W=1/N^2**

**> for i to N do**

**> for j to N do**

**> for k to N do**

**> for l to N do**

**> W[i,j,k,l] := 1/N^2;**

**> end do;**

**> end do;**

**> end do;**

**> end do;**

**> **

**> # Now compute the value for the derivative of B[i,j] from the B[i,j],**

**> # W[i,j,k,l], N, and t placing the values in the output array YP**

**> for i to N do**

**> for k to N do**

**> YP[N*(i-1)+k] := add(add(W[j,l,i,k]*B[j,l]**

**> -W[i,k,j,l]*B[i,k],l=1..N),j=1..N);**

**> end do;**

**> end do;**

**> end proc:**

3) Now at this point, you have a mapping of the variables, and a procedure

that can be used to compute the diff(B[i,j](t),t) for the i,j with the

W you program, so all we need now is initial conditions, and to perform

the call to generate the dsolve/numeric procedure.

**> ini := Array([seq(seq(i-j,j=1..N),i=1..N)]);**

**> dsn := dsolve(numeric,procedure=dproc,initial=ini,start=0,procvars=dvars);**

And the procedure can be used to compute the solutions:

**> dsn(0);**

**> dsn(1);**

**> **

This post generated using the online HTML conversion tool

Download the original worksheet