Question: How do I turn my numerical solution into a procedure and plot the resulting graph.

Hi,

I have a second order differential equation

d2y/dt2 = -6.478831125*sin(y)

to be solved numerically. I've successfully been able to solve it using the 4th order Runge-Kutta method, however it is not properly written as a procedure and I'm unsure of how to do this.

So far I have:

R:= -6.478831125
z[0]:= 0:
y[0]:= Pi/2:
h:= 0.01:
t:= 10:
for i from 0 to t-1 by 1 do
c0:= evalf(h*R*sin(y[i])):
k0:= evalf(h*(z[i])):
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]+(h/6)*(k0+(2*k1)+(2*k2)+k3)):
y[i+1]:= evalf(y[i]+(h/6)*(c0+(2*c1)+(2*c2)+c3)):
end do;

This outputs all values of c0,k0,c1... and z[1],y[1]... up to z[10],y[10] and I believe the solution to be correct. I'd like to make this into a procedure so that I can plot the 2d graph using DEplot but am unsure of how to do so. Please could somebody suggest a way this could be done?

As a side note, is it possible to print just the values of y[i+1] and z[i+1] without the c0,k0... values for each iteration?