Question: Cannot plot numerical ode

Can anyone help me?

> 

restart: with(plots):   

> 

h:=0.01: n:=6000: Digits:=18:

> 

t[0]:=0: x[0]:=0.2: y[0]:=0.1: epsilon:=20:

> 

F:=(x,y)->epsilon*(1-x^2)*y-x:

> 

begin:=time():

> 

for k from 0 to n do

> 

t[k+1]:=t[k]+h;

> 

x[k+1]:=x[k]+h*y[k];

> 

y[k+1]:=y[k]+h*F(x[k],y[k]);

> 

pt[k]:=[t[k],x[k],y[k]]:

> 

end do:

> 

cpu_time:=time()-begin;

> 

pointplot3d([seq(pt[j],j=0..n)],axes=normal,symbol=cross,symbolsize=8,color=red,labels=["t","x","y"],orientation=[-90,0]);

> 

sys:=diff(X(t),t)=Y(t),diff(Y(t),t)=F(X(t),Y(t));

> 

vars:={X(t),Y(t)}: ic:=X(0)=0.2,Y(0)=0.1:

> 

sol:=dsolve({sys,ic},vars,numeric,stepsize=h,method=classical[foreuler],output=listprocedure):

> 

odeplot(sol,[t,X(t)],0..h*n,axes=normal,style=line,numpoints=n,labels=["t","x"]);

> 

XX:=eval(X(t),sol): x1:=XX(1);  pt[100];

 

 

Download 09-1-1.mws

Please Wait...