#Sowmya Srinivasan hw13.txt 06 March 2014

#read C13.txt 

OtherJesusPi:=proc(n) local i,S,c,c1,j : 

c:=1 : 
S:=29 : 

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

S : 

end: 


####################################

RK4:=proc(f,x,y,x0,y0,x1,h) local xc,yc,k1,k2,k3,k4,L  : 

L:=[y0] : 

xc:=x0 : 
yc:=y0 : 

while xc<x1 do 
   k1:=subs({x=xc,y=yc},f) : 
   k2:=subs({x=xc+(1/2)*h,y=yc+(1/2)*h*k1},f) : 
   k3:=subs({x=xc+(1/2)*h,y=yc+(1/2)*h*k2},f) : 
   k4:=subs({x=xc+h,y=yc+h*k3},f) : 
   yc:=yc+(h/6)*(k1+2*k2+2*k3+k4) : 
   xc:=x0+h : 
   L:=[op(L),yc] : 
od : 

L : 

end: 


#########################################
#I was in a hurry so didn't get a chance to debug EulerH yet

EulerH:=proc(f,x,y,x0,y0,x1,Di) local i,h,err,u,v  : 

if DS(f,x,y,x0,y0)=FAIL then 
   RETURN(FAIL) : 
fi : 

i:=1: 
err:= 1 : 

while err>10^(-Di) do 
   h:=(x1-x0)/i : 
   u:=subs(x=x1,Euler(f,x,x0,y0,x1,h)[i+1] ) : 
   v:=subs(x=x1,DS(f,x,y,x0,y0)) : 
   err:=abs( u-v) : 
   i:=i+1 : 
od : 

h: 

end :