pds:-value(t=0.5)(0);
pds:-value(t=0.5)(1);

The above will give you values for t=0.5 with z=0 and then with z=1 and so on. 

Look at the help pages for the assign command:

?assign

For example: 

>ans:=solve( {x+y=1, 2*x+y=3}, {x, y} );

                    {x = 2, y = -1}

>assign(ans);

>x,y;

                     2, -1                   

restart:

f:=x/sqrt(1-cos(x));

limit(f,x=0,complex);

From the help files:

If dir is complex, the limit point infinity denotes complex infinity, that is, all infinities in the complex plane. If dir is real, the limit point infinity denotes both positive and negative infinity, and the limit is done bidirectionally. Otherwise, the limit point infinity denotes positive infinity, and -infinity denotes negative infinity.

restart: with(plots):

deq := (diff(y(x), x, x))*y(x)-1+x*y(x)^2 = cos(x);

sol := dsolve({deq, y(0) = 1, (D(y))(0) = 1}, numeric,range=-2..2);

p1:=odeplot(sol):

p2:=plot(x->op(2,op(2,sol(x))),-1..1,color=blue,style=point,symbol=circle):

display(p1,p2);

The plot looks like cos(x) over a small range.

This should clear it up a little.

plots:-display(DEplot(ode1, q(t), t = 0 .. 8760, q(t) = 0 .. 10000, [q(0) = 0], title = "contaminant accumulation and dissipation",linecolor=red,color=red), DEplot(ode2, q(t), t = 8760 .. 17520, q(t) = 0 .. 10000, [q(8760) = 9278],linecolor=blue,color=blue));

The method below works using odeplot.  

restart: with(plots):

p1:= dsolve({diff(q(t), t) = 3-3*q(t)/10^4, q(0)=0}, type=numeric, range=0..16000):

plot1:=odeplot(p1):

p2:=dsolve({diff(q(t), t) = -3*q(t)/10^4, q(0)=17520}, type=numeric, range=0..17520):

plot2:=odeplot(p2,color=blue):

display(plot1,plot2,title = "contaminant accumulation in the first year(red) and second year(blue)");
 

restart:

t:=[1,-1,0,1/2,3,3/2,sqrt(8)];

sort(t);

x:=evalf(t);

sort(x);

identify(%);

For your example as you posted it, you can do this after your t.

identify( sort(evalf(t)) );

Can you provide more details as to what you are trying to do, with an example, if possible?

restart:

x:=[Pi,sqrt(2), sqrt(3), sqrt(5), sqrt(10), sqrt(7)];

[x=Pi, 2^(1/2), 3^(1/2), 5^(1/2), 10^(1/2), 7^(1/2)]

identify(sort(evalf(x)));

[2^(1/2), 3^(1/2), 5^(1/2), 7^(1/2), Pi, 10^(1/2)]

This book looks like it covers the topics of intererest.  An advantage is the available code for you to try.

NonlinearFit from the Statistics package does not return a very good match to your data, using the model function y=ax^n+bx+c that you supplied.  Your might try another model function.  Otherwise, look into LSSolve, which is part of the Optimization package.

Download 221_Fit_from_data.mws

restart:

with(Student[Calculus1]):

Roots(cos(2*theta)=1/2,theta=0..2*Pi,numeric):

identify(%);

The help page for the Fit command has an example using X,Y data.

?Fit to bring it up.

Here is one way to do this without using the "sphere" command.

restart: with(plots):

animate( implicitplot3d, [x^2+y^2+z^2=A, x=-2..2, y=-2..2, z=-2..2], A=1..4 );

You can substitute your own equation for the above.  Look into the animate and the implicitplot3d commands for more details.  

Another way using rational values similar to what was done above in previous replies.

Download 221_solutions_lost.mws

The Maplesoft applications center has an example.  Here is the link.  The worksheet below is a modified version of the example from the link to make it look like the equation that you supplied.  The initial conditions are the ones from the example.

Download 221_vander.mws