10 Reputation

6 years, 163 days

Hi CarlI have attached my worksheet this...

Hi Carl

I have attached my worksheet this time, I am really struggling to see the error, the runge kutta formula only displays (tn,yn) I am having difficulty in obtaining the correct formula when including (x,y,z) I can not establish what I am doing different with each term ??

Thanks again for all your help this ones driving me crazy!

 > restart:
 >

Numerical Solution

 > with(plots):with(DEtools):
 >
 >

Equations & Initial conditions

 > de1:=D(x)(t)=a*x(t)-q*x(t)*y(t)-nn*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:=1.78:m:=0.00124:nn:=0.1:T:=10:N:=T/h:h:=0.01:
 > soln:=dsolve({de1,de2,de3,x(0)=38,y(0)=1,z(0)=2},{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..10,numpoints=100,title=Rabbit_v_fox_stoat):pexact:=%:
 > DEplot3d([de1,de2,de3],[x(t),y(t),z(t)],t=0..100,[[x(0)=38,y(0)=1,z(0)=2]],x=0..400,y=0..400,z=0..400,linecolor=black,arrows=medium,animatecurves=true,numframes=100):xyzexact:=%:
 (1)
 > soln(1);
 (2)

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-nn*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-nn*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]:=38:y[0]:=1:z[0]:=2:
 >

Step length and number of steps

 > h:=0.01:N:=T/h:
 >

Get the numerical solution

 > ESys(t,x,y,z,q,N):
 > meuler:=(x[100],y[100],z[100]);
 > 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=green):
 > ay:=plot({seq([t[n],y[n]],n=0..N)},style=point,colour=blue):
 > az:=plot({seq([t[n],z[n]],n=0..N)},style=point,colour=red):
 > axyz:=pointplot3d({seq([x[n],y[n],z[n]],n=0..N)},style=point,colour=black):
 >
 > display(array([display(pexact,ax,ay,az),display(xyzexact,axyz)])):
 (3)
 >

Initial conditions

 > t[0]:=0:x[0]:=38:y[0]:=1:z[0]:=2:
 >

Step length and number of steps

 > h:=0.01:N:=T/h:
 >

Get the numerical solution

 > Fourth(t,x,y,z,q,N):
 > four:=(x[100],y[100],z[100]);
 > 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=green):
 > iy:=plot({seq([t[n],y[n]],n=0..N)},style=point,colour=blue):
 > iz:=plot({seq([t[n],z[n]],n=0..N)},style=point,colour=red):
 > ixyz:=pointplot3d({seq([x[n],y[n],z[n]],n=0..N)},style=point,colour=black):
 >
 > display(array([display(pexact,ix,iy,iz),display(xyzexact,ixyz)])):
 >
 >
 (4)
 >
 >
 >
 >
 >

Hi Carl I have attached my worksheet thi...

Hi Carl

I have attached my worksheet this time, I am really struggling to see the error, the runge kutta formula only displays (tn,yn) I am having difficulty in obtaining the correct formula when including (x,y,z) I can not establish what I am doing different with each term ??

Thanks again for all your help this ones driving me crazy!

 > restart:
 >

Numerical Solution

 > with(plots):with(DEtools):
 >
 >

Equations & Initial conditions

 > de1:=D(x)(t)=a*x(t)-q*x(t)*y(t)-nn*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:=1.78:m:=0.00124:nn:=0.1:T:=10:N:=T/h:h:=0.01:
 > soln:=dsolve({de1,de2,de3,x(0)=38,y(0)=1,z(0)=2},{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..10,numpoints=100,title=Rabbit_v_fox_stoat):pexact:=%:
 > DEplot3d([de1,de2,de3],[x(t),y(t),z(t)],t=0..100,[[x(0)=38,y(0)=1,z(0)=2]],x=0..400,y=0..400,z=0..400,linecolor=black,arrows=medium,animatecurves=true,numframes=100):xyzexact:=%:
 (1)
 > soln(1);
 (2)

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-nn*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-nn*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]:=38:y[0]:=1:z[0]:=2:
 >

Step length and number of steps

 > h:=0.01:N:=T/h:
 >

Get the numerical solution

 > ESys(t,x,y,z,q,N):
 > meuler:=(x[100],y[100],z[100]);
 > 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=green):
 > ay:=plot({seq([t[n],y[n]],n=0..N)},style=point,colour=blue):
 > az:=plot({seq([t[n],z[n]],n=0..N)},style=point,colour=red):
 > axyz:=pointplot3d({seq([x[n],y[n],z[n]],n=0..N)},style=point,colour=black):
 >
 > display(array([display(pexact,ax,ay,az),display(xyzexact,axyz)])):
 (3)
 >

Initial conditions

 > t[0]:=0:x[0]:=38:y[0]:=1:z[0]:=2:
 >

Step length and number of steps

 > h:=0.01:N:=T/h:
 >

Get the numerical solution

 > Fourth(t,x,y,z,q,N):
 > four:=(x[100],y[100],z[100]);
 > 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=green):
 > iy:=plot({seq([t[n],y[n]],n=0..N)},style=point,colour=blue):
 > iz:=plot({seq([t[n],z[n]],n=0..N)},style=point,colour=red):
 > ixyz:=pointplot3d({seq([x[n],y[n],z[n]],n=0..N)},style=point,colour=black):
 >
 > display(array([display(pexact,ix,iy,iz),display(xyzexact,ixyz)])):
 >
 >
 (4)
 >
 >
 >
 >
 >

@Carl Love  Hi Carl, Many thanks fo...

Hi Carl,

Many thanks for the ongoing advice, not quite sure what happened with attachment I have pasted code below. In this code I have solved for t=1 with a step lenght of 0.01 hence x[100],y[100],z[100] for some unkown reason my fourth order estimate is further away than my mod euler. Is there something obvious to you that i have mistaken in my code.

Thanks again

Equations & Initial conditions

> de1:=D(x)(t)=a*x(t)-q*x(t)*y(t)-nn*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:=1.78:m:=0.00124:nn:=0.1:T:=10:N:=T/h:h:=0.01:

> soln:=dsolve({de1,de2,de3,x(0)=38,y(0)=1,z(0)=2},{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..10,numpoints=100,title=Rabbit_v_fox_stoat):pexact:=%:

DEplot3d([de1,de2,de3],[x(t),y(t),z(t)],t=0..100,[[x(0)=38,y(0)=1,z(0)=2]],x=0..400,y=0..400,z=0..400,linecolor=black,arrows=medium,animatecurves=true,numframes=100):xyzexact:=%:

> soln(1);                                                                  t=1 x=2829.68  y=0.61  z=0.75

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-nn*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-nn*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]:=38:y[0]:=1:z[0]:=2:

Step length and number of steps

> h:=0.01:N:=T/h:

Get the numerical solution

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

> (x[100],y[100],z[100]);                                                                x=2587.22  y=0.59  z=0.69

> 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=green):

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

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

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

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

Initial conditions

> t[0]:=0:x[0]:=38:y[0]:=1:z[0]:=2:

>

Step length and number of steps

> h:=0.01:N:=T/h:

>

Get the numerical solution

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

> (x[100],y[100],z[100]);                                                             x=1206.79   y=0.52    z=0.52

> 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=green):

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

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

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

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

@Carl Love  Hi Carl, Many thanks fo...

Hi Carl,

Many thanks for the ongoing advice, not quite sure what happened with attachment I have pasted code below. In this code I have solved for t=1 with a step lenght of 0.01 hence x[100],y[100],z[100] for some unkown reason my fourth order estimate is further away than my mod euler. Is there something obvious to you that i have mistaken in my code.

Thanks again

Equations & Initial conditions

> de1:=D(x)(t)=a*x(t)-q*x(t)*y(t)-nn*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:=1.78:m:=0.00124:nn:=0.1:T:=10:N:=T/h:h:=0.01:

> soln:=dsolve({de1,de2,de3,x(0)=38,y(0)=1,z(0)=2},{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..10,numpoints=100,title=Rabbit_v_fox_stoat):pexact:=%:

DEplot3d([de1,de2,de3],[x(t),y(t),z(t)],t=0..100,[[x(0)=38,y(0)=1,z(0)=2]],x=0..400,y=0..400,z=0..400,linecolor=black,arrows=medium,animatecurves=true,numframes=100):xyzexact:=%:

> soln(1);                                                                  t=1 x=2829.68  y=0.61  z=0.75

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-nn*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-nn*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]:=38:y[0]:=1:z[0]:=2:

Step length and number of steps

> h:=0.01:N:=T/h:

Get the numerical solution

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

> (x[100],y[100],z[100]);                                                                x=2587.22  y=0.59  z=0.69

> 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=green):

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

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

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

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

Initial conditions

> t[0]:=0:x[0]:=38:y[0]:=1:z[0]:=2:

>

Step length and number of steps

> h:=0.01:N:=T/h:

>

Get the numerical solution

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

> (x[100],y[100],z[100]);                                                             x=1206.79   y=0.52    z=0.52

> 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=green):

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

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

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

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

pointplot3d and global/local conflict...

Hi Carl

Many thanks for the adivce the code and plot seem to work brilliantly now! I have one another small question for you too.

In my code I am looking to obtain a value for a given (t) for my differential equations and also my Esys and Fourth estimates. I obtained the values for t=5 by using the following command;

soln(5);

how could I extract thesesvalues form my estimates ( Esys and Fourth ) ?

Best Regards

Anthony

ps (file attached)

pointplot3d and global/local conflict...

Hi Carl

Many thanks for the adivce the code and plot seem to work brilliantly now! I have one another small question for you too.

In my code I am looking to obtain a value for a given (t) for my differential equations and also my Esys and Fourth estimates. I obtained the values for t=5 by using the following command;

soln(5);

how could I extract thesesvalues form my estimates ( Esys and Fourth ) ?

Best Regards

Anthony

ps (file attached)

 Page 1 of 1
﻿