#hw13.txt, 3/9/14, Anthony Zaleski

Help:=proc():
print(`OtherJesusPi(n),RK4(f,x,y,x0,y0,x1,h)`):
end:


##########Problem 1: Jesus' second formula##########
#The nth partial sum given in the second formula
#in https://sites.google.com/site/guilleramath/

OtherJesusPi:=proc(n) local S,i,j,a,b:

S:=0:
a:=1: b:=2880^3:

for i from 0 to n do
	S:=S+a*(5418*i^2+693*i+29):
	a:=-a*mul(6*i+j,j=1..6)/b/(i+1)^6:
od:

S:
end:

#Output:

#evalf(128*sqrt(5)/Pi^2,25);
#28.99981493568089545672014

#for i from 1 to 4 do evalf(OtherJesusPi(1),25); od;
#28.99981493537808641975309
#28.99981493568089596882878
#28.99981493568089545671924
#28.99981493568089545672013

##########Problem 2: Runge-Kutta##########

#One step of R-K:
RK1:=proc(f,x,y,x0,y0,h) local k1,k2,k3,k4:
k1:=subs({x=x0,y=y0},f):
k2:=subs({x=x0+h/2,y=y0+h*k1/2},f):
k3:=subs({x=x0+h/2,y=y0+h*k2/2},f):
k4:=subs({x=x0+h/2,y=y0+h*k2/2},f):
k4:=subs({x=x0+h,y0+h*k3},f):

y0+h*(k1+2*k2+2*k3+k4)/6:
end:



#RK4(f,x,y,x0,y0,x1,h): R-K method approximating
#y'(x)=f(x,y), y(x0)=y0, all the way to x1
#with increments of h


RK4:=proc(f,x,y,x0,y0, x1, h) local L,xc,yc:
L:=[y0]:

xc:=x0:
yc:=y0:

while xc<x1 do

yc:=RK1(f,x,y,xc,yc,h):
xc:=xc+h:
L:=[op(L),yc]:

od:
L:

end: