Fantastic, that's just what I was looking for! I figured out shortly after posting that I could use pointplot3d to plot the individual points but couldn't figure out how to plot the line between points. Thanks Carl.

Fantastic, that's just what I was looking for! I figured out shortly after posting that I could use pointplot3d to plot the individual points but couldn't figure out how to plot the line between points. Thanks Carl.

In regards to what I posted earlier I think I have realised one of my mistakes. For a step size of 0.01 and a period of 10 seconds I believe that I should have actually used N=1000.

As well as the phase plane diagram, I am trying to plot the values of z and y against time on the same graph. For this I introduced t as seen below however this produces a long error when I attempt to plot just y against t.

restart;
Digits:=15:
R:= -6.478831125:
h:= 0.01:
N:= 1000:
z:=Vector(N+1):
y:=Vector(N+1):
t:=Vector(N+1):
z[1]:= 0:
y[1]:= Pi/2:
t[1]:= 0:
for i from 1 to N do
c0:= evalf(h*R*sin(y[i])): #increment for z
k0:= evalf(h*(z[i])): #increment for y
c1:= evalf(h*R*sin(y[i]+(k0/2))):
k1:= evalf(h*(z[i]+(c0/2))):
c2:= evalf(h*R*sin(y[i]+(k1/2))):
k2:= evalf(h*(z[i]+(c1/2))):
c3:= evalf(h*R*sin(y[i]+k2)):
k3:= evalf(h*(z[i]+c2)):
z[i+1]:= evalf(z[i]+(1/6)*(c0+(2*c1)+(2*c2)+c3)):
y[i+1]:= evalf(y[i]+(1/6)*(k0+(2*k1)+(2*k2)+k3)):
t[i+1]:= evalf(t[i]+h):
end do:
plot(y,z);
plot(t,Y);

In regards to what I posted earlier I think I have realised one of my mistakes. For a step size of 0.01 and a period of 10 seconds I believe that I should have actually used N=1000.

As well as the phase plane diagram, I am trying to plot the values of z and y against time on the same graph. For this I introduced t as seen below however this produces a long error when I attempt to plot just y against t.

restart;
Digits:=15:
R:= -6.478831125:
h:= 0.01:
N:= 1000:
z:=Vector(N+1):
y:=Vector(N+1):
t:=Vector(N+1):
z[1]:= 0:
y[1]:= Pi/2:
t[1]:= 0:
for i from 1 to N do
c0:= evalf(h*R*sin(y[i])): #increment for z
k0:= evalf(h*(z[i])): #increment for y
c1:= evalf(h*R*sin(y[i]+(k0/2))):
k1:= evalf(h*(z[i]+(c0/2))):
c2:= evalf(h*R*sin(y[i]+(k1/2))):
k2:= evalf(h*(z[i]+(c1/2))):
c3:= evalf(h*R*sin(y[i]+k2)):
k3:= evalf(h*(z[i]+c2)):
z[i+1]:= evalf(z[i]+(1/6)*(c0+(2*c1)+(2*c2)+c3)):
y[i+1]:= evalf(y[i]+(1/6)*(k0+(2*k1)+(2*k2)+k3)):
t[i+1]:= evalf(t[i]+h):
end do:
plot(y,z);
plot(t,Y);

Thank you very much for your reply.

Could you explain why we should define y and z as vectors?

Edit: I have realised some of my errors, reply below.

Thank you very much for your reply.

Could you explain why we should define y and z as vectors?

Edit: I have realised some of my errors, reply below.

Fantastic, thanks very much for the help Carl :)

Fantastic, thanks very much for the help Carl :)

Thanks for your reply Carl, I really appreciate it.

After reading your comment I've taken another look at my procedures and made some changes. For Simpson's rule I think I may have found 4 of the errors:

> f:= x -> cos(x)*exp(-x/4):

> CompTrap:= proc(f,a,b,n)
local i, aa, h;
h:= (b-a)/n;
aa:= f(a)+f(b);

for i from 1 to n-1 by 1 do

aa:= aa + 2*(f(a+(i*h)));
end do;
aa:= aa*(1/2)*h;
return aa;
end proc:

> Comp_Simp_1_3:= proc(f,a,b,n)
local i, aa, h;
h:= (b-a)/n;
aa:= f(a)+f(b);

for i from 1 to n-1 by 1 do

aa:= aa + 4*f(a+(i-1)*h) + 2*f(a+(i*h));

end do;
aa:= aa*(1/3)*h;
return aa;
end proc:


>evalf(int(cos(x)*(exp(-x/4)), x=1..8));
-.3871279027

> evalf(CompTrap(f,1,8,8));
-.3462142474

> evalf(Comp_Simp_1_3(f,1,8,8));
-.5663721786

My CompTrap now converges to the calculated value when n is increased, however CompSimp is converging towards a different value so there are still mistakes in my procedure that I'm having difficulty spotting.

Thanks for your reply Carl, I really appreciate it.

After reading your comment I've taken another look at my procedures and made some changes. For Simpson's rule I think I may have found 4 of the errors:

> f:= x -> cos(x)*exp(-x/4):

> CompTrap:= proc(f,a,b,n)
local i, aa, h;
h:= (b-a)/n;
aa:= f(a)+f(b);

for i from 1 to n-1 by 1 do

aa:= aa + 2*(f(a+(i*h)));
end do;
aa:= aa*(1/2)*h;
return aa;
end proc:

> Comp_Simp_1_3:= proc(f,a,b,n)
local i, aa, h;
h:= (b-a)/n;
aa:= f(a)+f(b);

for i from 1 to n-1 by 1 do

aa:= aa + 4*f(a+(i-1)*h) + 2*f(a+(i*h));

end do;
aa:= aa*(1/3)*h;
return aa;
end proc:


>evalf(int(cos(x)*(exp(-x/4)), x=1..8));
-.3871279027

> evalf(CompTrap(f,1,8,8));
-.3462142474

> evalf(Comp_Simp_1_3(f,1,8,8));
-.5663721786

My CompTrap now converges to the calculated value when n is increased, however CompSimp is converging towards a different value so there are still mistakes in my procedure that I'm having difficulty spotting.