It is not necessary to store all the (possibly many and large) procedures generated by instantiating at different numeric values for parameter **delta**, in order to get decent speed-up by avoiding setting the parameter **delta** inefficiently.

That requires later clean-up (**forget**) for the memory to be cleared, and for some examples is an extra potential allocation issue. Such memory concerns and the use of advanced **printf** functionality are unnecessary for accomplishing this plotting task reasonably efficiently. The OP has made no request to store such procedures externally, and that's a pretty rare task.

Instead of the above, either name (here **t** or **delta**) can be put on the first axis in the 3D plot, while still getting the speedup of setting values efficiently for parameter **delta**. The *parametric calling sequence* of the **plot3d** command allows for parameter **delta** to be used more efficiently as the first ("outer") range, while still flexibly allowing either **t,delta** or **delta,t** being ascribed to axis[1],axis[2] respectively.

This below is done in Maple 2017.3, since the OP was using that version a while back.

restart;
phi:=0:lambda:=0.1:N:=5:M:=sqrt(N*(N+1))*exp(I*phi):omegap:=10:
var:={n(t),u(t),z(t)}:
dsys:={diff(n(t),t)=-2*(n(t)-N)+(u(t)-M)*exp(-2*I*omegap*t/lambda)
+((z(t)-conjugate(M))*exp(2*I*omegap*t/lambda)),
diff(u(t),t)=-2*(1-I*delta)*u(t)
+2*(n(t)-N)*exp(2*I*omegap*t/lambda)+2*M,
diff((z(t),t))=-2*(1+I*delta)*z(t)
+2*(n(t)-N)*exp(-2*I*omegap*t/lambda)
+2*conjugate(M)}:
res1:=dsolve(dsys union {n(0)=0,u(0)=0,z(0)=0},
numeric, parameters=[delta], output=listprocedure):
Nt:=eval(n(t),res1):
Ntd:=proc(t,delta)
if not [args]::list(numeric) then 'procname'(args);
else
if rhs~(Nt(parameters))<>[delta] then
Nt(parameters=[delta]);
end if;
Nt(t);
end if;
end proc:
CodeTools:-Usage(plot3d([t, delta, Ntd(t,delta)],
delta=-10..10, t=0..1, axes=boxed));
memory used=8.44GiB, alloc change=36.00MiB, cpu time=64.89s,
real time=64.91s, gc time=5.06s

CodeTools:-Usage(plot(Ntd(.5,delta), delta=-10..10,
numpoints=49, adaptive=false));
memory used=4.66GiB, alloc change=80.00MiB, cpu time=35.69s,
real time=34.11s, gc time=3.65s

[edit] I wrote the above for Maple 2017.2, but in Maple 2019.2 (and later) some improvement to **dsolve**'s runtime parameter handling allows for even simpler and more understandable code. The conditional test to check whether the parameter value has changed is no longer needed above in order to get the speedup from that use of the parametric call to **plot3d**. So the following is enough, paring it down while keeping it understandable (IMNSHO):

phi:=0:lambda:=0.1:N:=5:M:=sqrt(N*(N+1))*exp(I*phi):omegap:=10:
dsys:={diff(n(t),t)=-2*(n(t)-N)+(u(t)-M)*exp(-2*I*omegap*t/lambda)
+((z(t)-conjugate(M))*exp(2*I*omegap*t/lambda)),
diff(u(t),t)=-2*(1-I*delta)*u(t)
+2*(n(t)-N)*exp(2*I*omegap*t/lambda)+2*M,
diff((z(t),t))=-2*(1+I*delta)*z(t)
+2*(n(t)-N)*exp(-2*I*omegap*t/lambda)
+2*conjugate(M)}:
res1:=dsolve(dsys union {n(0)=0,u(0)=0,z(0)=0},
numeric, parameters=[delta], output=listprocedure):
Nt:=eval(n(t),res1):
Ntd:=proc(t,delta) Nt(parameters=[delta]); Nt(t); end proc:
plot3d([t, delta, 'Ntd'(t,delta)], delta=-10..10, t=0..1);