The code below is nearly five times faster than the one I posted yesterday (see above).
It uses only slightly more memory.

The tip from Preben - to convert from differential operators to 'D' and then substitute to an inert name - saved an unnecessary loop and could be of value for many Maple users, dealing with solution of equations containing derivatives.

Now the code is fast enough.
So - thank you Preben, for all the good advice.

If anyone can figure out why I need to 'clean' f by writing it and subsequently reading it from a file, I would be most interested ;) 

restart;
#printlevel:=3:
Ns:=16:
NP:=50:
NC:=10:
NS:=NP+NC:

Z:=Array(1..Ns,1..NS);
r:=rand(0..99);
r();

for s from 1 to Ns do
for i from 1 to NS do
Z[s,i]:=1.0*r();
od: od:

for s from 1 to Ns do
for i from 1 to NP do
phi[s][i]:=add(z[s][j]^2,j=1..NS)
od: od;

for s from 1 to Ns do
for i from 1 to NP do
f[s][i]:=expand(z[s][i]-phi[s][i]);
od: od:

for s from 1 to Ns do
alias(seq(z[s][i]=z[s][i](seq(z[s][j],j=NP+1..NP+NC)),i=1..NP));
od:

save f,"File":
unassign('f'):
read "File":

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
df_dz[s][i,j]:=diff(f[s][i],z[s][j]);
df_dz[s][i,j]:=subs(D=DD,convert(df_dz[s][i,j],D));
od: od: od:

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
DD[j-NP](z[s][i])(seq(z[s][k],k=NP+1..NP+NC)):=Dpc[s][i,j];
od: od: od:

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
df_dz[s][i,j]:=eval(df_dz[s][i,j]);
od: od: od:

for s from 1 to Ns do
for i from 1 to NS do
z[s][i]:=Z[s][i];
od: od:

for s from 1 to Ns do
for j from NP+1 to NP+NC do
Sol:=solve({seq(df_dz[s][i,j],i=1..NP)},{seq(Dpc[s][k,j],k=1..NP)});
assign(Sol):
od: od:

Unfortunately the code below, which is not very sophisticated, runs faster and uses significantly less memory.
This holds, in particular, as Ns is increased; I have used the cases (NP,NC)=(50,10) and Ns=2,4 and 8.

My original problem with the code was that nearly all CPU time is used for the loop with
df_dz[s][i,j]:=eval(df_dz[s][i,j],diff(z[s][p](seq(z[s][j],j=NP+1..NP+NC)),z[s][j])=Dpc[s][p,j]);
I turned to this forum for advice, since I believed that there should exist a more efficient way.

The Dpc variable is the one I finally want to solve for (it is used elsewhere in the main code), but I introduced it rather clumsily in this loop.

The reason I built this loop was that I could not assign Z values to z in the expressions for df_dz at this stage without running into problems with the diffferential operators. (Z values must be given before 'solve' is applied; they are computed elsewhere in the main code but are in this example obtained from a random number generator).

I have tried to follow the advice given here to switch to the 'D' operator somehow so that Z values could be assigned, but my Maple skills were too limited; I could not get this working.

I have a vague feeling that the code below could be made substantially more efficient. Maple syntax is, however, not always intuitive and I cannot presently see how this could be done. 

restart;
#printlevel:=3:
Ns:=8:
NP:=50:
NC:=10:
NS:=NP+NC:

Z:=Array(1..Ns,1..NS);
r:=rand(0..99);
r();

for s from 1 to Ns do
for i from 1 to NS do
Z[s,i]:=1.0*r();
od: od:

for s from 1 to Ns do
for i from 1 to NP do
phi[s][i]:=add(z[s][j]^2,j=1..NS)
od: od;

for s from 1 to Ns do
for i from 1 to NP do
f[s][i]:=expand(z[s][i]-phi[s][i]);
od: od:

for s from 1 to Ns do
alias(seq(z[s][i]=z[s][i](seq(z[s][j],j=NP+1..NP+NC)),i=1..NP));
od:

save f,"File":
unassign('f'):
read "File":

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
df_dz[s][i,j]:=diff(f[s][i],z[s][j]);
od: od: od:

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
for p from 1 to NP do
df_dz[s][i,j]:=eval(df_dz[s][i,j],diff(z[s][p](seq(z[s][j],j=NP+1..NP+NC)),z[s][j])=Dpc[s][p,j]);
od:
od: od: od:

for s from 1 to Ns do
for i from 1 to NS do
z[s][i]:=Z[s][i];
od: od: 

for s from 1 to Ns do
for j from NP+1 to NP+NC do
Sol:=solve({seq(df_dz[s][i,j],i=1..NP)},{seq(Dpc[s][k,j],k=1..NP)});
assign(Sol):
od: od:

Thanks again for your help!

However, in my code Z is a real variable, but for real values of Z the solution you suggested unfortunately does not work.

In the code below, for example, the assignment Z[s,i]:=1*r() does work, whereas Z[s,i]:=1.0*r() does not. 
This really puzzles me; how can this be the case?

You may also note that my implementation below for assigning values to Dpc is different than the one you suggested in your last reply. This is because it is much more efficient for large values of Ns, NP and NC.
The code below runs a factor 7 faster (with slightly less memory use) for (Ns,NP,NC)=(5,50,10).
Its CPU time is nearly linear with Ns, whereas the CPU time of your solution is rather strongly dependent on Ns.

Again, these runs were made for Z[s,i]:=1*r(); since real Z values fail for both implementations.

Any help here would be really appreciated :) 

restart;
printlevel:=3:
Ns:=2:
NP:=2: #=(m_max+1)*(n_max-1) = 6*8 = 50
NC:=2: #NC/NP=2/(n_max-1) = 2/10 = 1/5
NS:=NP+NC:

Z:=Array(1..Ns,1..NS);
#Z:=Array(1..Ns,1..NS,(s,i)->s*i);
r:=rand(0..99);
r();

for s from 1 to Ns do
for i from 1 to NS do
Z[s,i]:=1*r();
od: od:

for s from 1 to Ns do
indepvar[s]:=seq(z[s][j],j=NP+1..NP+NC);
depvar[s]:=seq(z[s][j](indepvar[s]),j=1..NP);
var[s]:=depvar[s],indepvar[s]
end do;

for s from 1 to Ns do
for i from 1 to NP do
phi[s][i]:=add(var[s][j]^2,j=1..NS)
od: od;

for s from 1 to Ns do
for i from 1 to NP do
f[s][i]:=expand(depvar[s][i]-phi[s][i]);
od: od;

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
df_dz[s][i,j]:=diff(f[s][i],z[s][j]);
df_dz[s][i,j]:=convert(df_dz[s][i,j],D)
od: od: od;

vals:={seq(seq(var[s][i]=Z[s,i],i=1..NS),s=1..Ns)};
eqs:=eval({seq(seq(seq(df_dz[s][i,j],i=1..NP),j=NP+1..NP+NC),s=1..Ns)},vals):
indets(eqs,function):
Sol:=solve(eqs,indets(eqs,function)):

assign(Sol):
assign(vals):

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
Dpc[s][i,j]:=D[j-NP](z[s][i])(seq(z[s][k],k=NP+1..NP+NC));
#print(s,i,j,Dpc[s][i,j]):
od: od: od;

Thanks again, Preben!
You really find elegant syntax solutions; I am impressed.

For me, the most important thing you have shown is how to substitute values into formulas containing derivatives; you do this by using convert(...,D). 
You have now also shown how variable dependencies can be expressed without complicating matters using 'alias'.

As seen above, I want the solution to the equations to be stored in the array Dpc.
Would the coding below be an efficient way to store the solution?

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
Dpc[s][i,j]:=D[j-NP](z[s][i])(seq(z[s][k],k=NP+1..NP+NC));
od: od: od;

Finally, I have one remaining problem: the following line

vals:={seq(seq(var[s][i]=s*i,i=1..NS),s=1..Ns)};

works well for this simplified test case, but when I, as in my code, want to use previously stored values from an array;

vals:={seq(seq(var[s][i]=Z[s][i],i=1..NS),s=1..Ns)};

solve fails.
What can possibly be the reason for this?
Any help would be gratefully received :)

Thanks a lot, Preben, your swift and knowledgeable answer has been really helpful.
As you may have guessed, my question was a simplification of the actual Maple coding problem that I have.

It is, after some irrelevant simplifications, more like the following.
Unfortunately these code lines are really inefficient and take a lot of CPU time.
My goal is to solve the resulting system of equations for the derivatives Dpc[s][k,j].

You can see that I use some unconventional paths, like writing on an external file, in order to get the syntax working.
Clearly there should be smoother and more efficient ways to write this piece of code.
I guess, as you suggest, that it is a better idea to use 'unapply' than 'alias' for a simpler syntax.
If I could get additional help with simplifications here, I would be extremely grateful:)

restart:
Ns:=1:
NP:=2:
NC:=2:
NS:=NP+NC:

for s from 1 to Ns do
for i from 1 to NP do
phi[s][i]:=sum(z[s][j]**2,j=1..NS);
od: od:

for s from 1 to Ns do
for i from 1 to NP do
f[s][i]:=expand(z[s][i]-phi[s][i]);
od: od:

for s from 1 to Ns do
alias(seq(z[s][i]=z[s][i](seq(z[s][j],j=NP+1..NP+NC)),i=1..NP));
od:

save f,"File":
unassign('f'):
read "File":

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
df_dz[s][i,j]:=diff(f[s][i],z[s][j]);
od: od: od:

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
for p from 1 to NP do
df_dz[s][i,j]:=eval(df_dz[s][i,j],diff(z[s][p](seq(z[s][j],j=NP+1..NP+NC)),z[s][j])=Dpc[s][p,j]); od: #print(s,i,j,df_dz[s][i,j]):
od: od: od:

for s from 1 to Ns do
for i from 1 to NS do
z[s][i]:=s*i;
od: od:

for s from 1 to Ns do
for j from NP+1 to NP+NC do
Sol:=solve({seq(df_dz[s][i,j],i=1..NP)},{seq(Dpc[s][k,j],k=1..NP)});
assign(Sol):
od: od:

for s from 1 to Ns do
for i from 1 to NP do
for j from NP+1 to NP+NC do
print(s,i,j,Dpc[s][i,j]);
od: od: od: