solution of system odes



Dear Users of the Forum,

I have a problem with a program which I would like to solve system of
odes, originating form method of lines. There is my program:

kill(all);
N:2;h:1/N; zz:2; fun(x):=1;
for i:1 while i<N+1 do
(
   eq:'diff(u,x)+u*(N+zz)-N*fun(x)=0,
   kill(ro),
   kill(rs),
   kill(zz1),
   gs:ode2(eq,u,x),
   ps:ic1(gs,x=0,u=exp(-i*h)),
   zz1:(zz-h*(ev(ps,x=2)-exp(-i*h)))/(1+h*zz),
   zz:rhs(float(zz1)),
   fun(x):=rhs(ps),
   plot2d(fun(x),[x,0,2]),
   kill(eq)
);

When N=2, then I get two plots. When N=3, then it doesn't produce the
last plot. I suppose that Maxima tends to solve the differential
equations exactly and it stops (or works extremely long) due to
complicated formula in 3rd equation. Is it possible to made Maxima
to use float point representation of terms such as %e^{...}, which aper
in solutions of odes? If finally my program works correctly, would it be
possible to gather all plots of solutions of odes at one 3d plot? I
would be grateful for any comments and hints.
Piotr.