Question: output of numerical solutions of ODE as a table

Hello experts in Maple,

I'm going to solve IVP numerically by this command bellow:

var:={u(t),v(t),w(t)}: 
dsys1 :=diff(u(t),t)=-(Gamma+AM*cos(q*t+phi))*u(t)+(Delta-AM*sin(q*t+phi))*v(t)+Omega*w(t), diff(v(t),t)=-(Delta+AM*sin(q*t+phi))*u(t)-(Gamma-AM*cos(q*t+phi))*v(t), diff(w(t),t)=-0.5-2*Gamma*w(t)-Omega*u(t):
dsol1 :=dsolve({dsys1,ini1},numeric,method=lsode,var, abserr=1e-9, relerr=1e-8,output=listprocedure):
dsolu:=subs(dsol1,u(t)):dsolv:=subs(dsol1,v(t)):dsolw:=subs(dsol1,w(t)):
t1:=array(0..N1,[]): u1:=array(0..N1,[]): v1:=array(0..N1,[]): w1:=array(0..N1,[]): pt1:=array(0..N1,[]):pt2:=array(0..N1,[]):pt3:=array(0..N1,[]): 
for i from 0 to N1 do t1[i]:=evalf(th*i):u1[i]:=evalf(dsolu(t1[i]));v1[i]:=evalf(dsolv(t1[i])):w1[i]:=evalf(dsolw(t1[i])):pt1[i]:=[t1[i],u1[i]]:pt2[i]:=[t1[i],v1[i]]:pt3[i]:=[t1[i],w1[i]]: od:
> 
unassign('i'):mytab1:=[seq(pt1[i],i=0..N1)]:mytab2:=[seq(pt2[i],i=0..N1)]:mytab3:=[seq(pt3[i],i=0..N1)]:

Now I would like to get table output. I know I shoud use fopen- fprintf-fclose.. in the loop but I did not know how??

please need advise.

Thankx in advance.