Question: could anyone help me with my euler plot too!

[I split this off as a separate Question, cleaned up the machine-generated Title, added the Tags, and cleaned out some excessive blank lines.---Carl Love as Moderator]

could anyone help me with my euler plot too!

> with(plots):with(DEtools):

Equations & Initial conditions

> de1:=D(x)(t)=a*x(t)-q*x(t)*y(t)-n*x(t)*z(t);

> de2:=D(y)(t)=-b*y(t)+k*x(t)*y(t);

> de3:=D(z)(t)=-c*z(t)+m*x(t)*z(t);

>

> a:=4.43:q:=0.022:b:=0.89:k:=0.00062:c:=0.25:m:=0.25:n:=1:T:=4:h:=0.2:N:=T/h:

> soln:=dsolve({de1,de2,de3,x(0)=10,y(0)=10,z(0)=10},{x(t),y(t),z(t)},numeric):

> odeplot(soln,[[t,x(t),color=green],[t,y(t),color=blue],[t,z(t),color=red]],0..100,numpoints=100,title=Rabbit_v_fox_stoat):pexact:=%:

> DEplot3d([de1,de2,de3],[x(t),y(t),z(t)],t=0..100,[[x(0)=10,y(0)=10,z(0)=10]],x=0..20,y=0..20,z=0..20,linecolor=black,arrows=medium,animatecurves=true,numframes=100):xyzexact:=%:

My numerical solve (Euler's method & 4th order)

> ESys:=proc(t,x,y,z,q,N)

> local f,g,e,n,k1,l1,m1;

> f:=(x,y,z)->a*x-q*x*y-n*x*z;

> g:=(x,y,z)->-b*y+k*x*y;

> e:=(x,y,z)->-c*z+m*x*z;

> for n from 0 to N do

> k1:=h*f(x[n],y[n],z[n]);

> l1:=h*g(x[n],y[n],z[n]);

> m1:=h*e(x[n],y[n],z[n]);

> t[n+1]:=t[n]+h;

> x[n+1]:=x[n]+k1;

> y[n+1]:=y[n]+l1;

> z[n+1]:=z[n]+m1;

> od;

> end proc:

>

> Fourth:=proc(t,x,y,z,q,N)

> local f,g,e,n,k1,k2,k3,k4,l2,l1,l3,l4,m1,m2,m3,m4;

> f:=(x,y,z)->a*x-q*x*y-n*x*z;

> g:=(x,y,z)->-b*y+k*x*y;

> e:=(x,y,z)->-c*z+m*x*z;

> for n from 0 to N do

> k1:=h*f(x[n],y[n],z[n]);

> l1:=h*g(x[n],y[n],z[n]);

> m1:=h*e(x[n],y[n],z[n]);

> k2:=h*f(x[n]+0.5*h,y[n]+0.5*k1,z[n]+0.5*k1);

> l2:=h*g(x[n]+0.5*h,y[n]+0.5*l1,z[n]+0.5*l1);

> m2:=h*e(x[n]+0.5*h,y[n]+0.5*m1,z[n]+0.5*m1);

> k3:=h*f(x[n]+0.5*h,y[n]+0.5*k2,z[n]+0.5*k2);

> l3:=h*g(x[n]+0.5*h,y[n]+0.5*l2,z[n]+0.5*l2);

> m3:=h*e(x[n]+0.5*h,y[n]+0.5*m2,z[n]+0.5*m2);

> k4:=h*f(x[n]+h,y[n]+k3,z[n]+k3);

> l4:=h*g(x[n]+h,y[n]+l3.z[n]+l3);

> m4:=h*e(x[n]+h,y[n]+m3,z[n]+m3);

> t[n+1]:=t[n]+h;

> x[n+1]:=x[n]+((1)/(6))*(k1+2*k2+2*k3+k4);

> y[n+1]:=y[n]+((1)/(6))*(l1+2*l2+2*l3+l4);

> z[n+1]:=y[n]+((1)/(6))*(m1+2*m2+2*m3+m4);

> od;

> end proc:

Initial conditions

> t[0]:=0:x[0]:=10:y[0]:=10:z[0]:=10:

Step length and number of steps

> h:=0.2:N:=T/h:

Get the numerical solution

> ESys(t,x,y,z,q,N):

> seq([t[n],evalf(x[n]),evalf(y[n]),evalf(z[n])],n=0..N):

Plot against the Maple's numerical solution

> ax:=plot({seq([t[n],x[n]],n=0..N)},style=point,colour=black):

> ay:=plot({seq([t[n],y[n]],n=0..N)},style=point,colour=black):

> az:=plot({seq([t[n],z[n]],n=0..N)},style=point,colour=black):

> axyz:=plot({seq([x[n],y[n],z[n]],n=0..N)},style=point,colour=black):

>

> display(array([display(pexact,ax,ay,az),display(xyzexact,axyz)]));

Error, (in plot) found points with fewer or more than 2 components

Error, (in plots:-display) expecting plot structure but received: axyz

Initial conditions

> t[0]:=0:x[0]:=10:y[0]:=10:z[0]:=10:

Step length and number of steps

> h:=0.2:N:=T/h:

Get the numerical solution

> Fourth(t,x,y,z,q,N):

> seq([t[n],evalf(x[n]),evalf(y[n]),evalf(z[n])],n=0..N):

Plot against the Maple's numerical solution

> ix:=plot({seq([t[n],x[n]],n=0..N)},style=point,colour=black):

> iy:=plot({seq([t[n],y[n]],n=0..N)},style=point,colour=black):

> iz:=plot({seq([t[n],z[n]],n=0..N)},style=point,colour=black):

> ixyz:=plot({seq([x[n],y[n],z[n]],n=0..N)},style=point,colour=black):

>

> display(array([display(pexact,ix,iy,iz),display(xyzexact,ixyz)]));

Error, (in plot) found points with fewer or more than 2 components

Error, (in plots:-display) expecting plot structure but received: ixyz