antfaulkner

10 Reputation

One Badge

5 years, 345 days

MaplePrimes Activity


These are replies submitted by antfaulkner

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:=%:

(D(x))(t) = a*x(t)-q*x(t)*y(t)-nn*x(t)*z(t)

 

(D(y))(t) = -b*y(t)+k*x(t)*y(t)

 

(D(z))(t) = -c*z(t)+m*x(t)*z(t)

(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)

 

 

 

 

 

 

Download NUMERICAL_MASTER_-_C.mw

@Carl Love 

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:=%:

(D(x))(t) = a*x(t)-q*x(t)*y(t)-nn*x(t)*z(t)

 

(D(y))(t) = -b*y(t)+k*x(t)*y(t)

 

(D(z))(t) = -c*z(t)+m*x(t)*z(t)

(1)

soln(1);

[t = 1., x(t) = HFloat(2829.675175405923), y(t) = HFloat(0.6118299784863152), z(t) = HFloat(0.7486718921094427)]

(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)])):

2587.223195, .5895973518, .6921971223

(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)])):

 

 

1206.792883, .5232528524, .5224305498

(4)

 

 

 

 

 


Download NUMERICAL_MASTER_-_C.mw

@Carl Love 

@Carl Love 

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 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)])):

 

 

 

 

 

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)

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