#Sowmya Srinivasan hw18.txt 31 March 2014 

#read C18.txt 


#SIrecS(L,Ini,f,m,N) outputs a list of N+1 SIrec(L,Ini,f,m,n) 

SIrecS:=proc(L,Ini,f,m,N) local n:
[ seq(SIrec(L,Ini,f,m,n),n=0..N)]:
end:


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


#dBVPIc(L,Ini,Fini,f,m,N) cleverly solves boundary value recurrence equations
dBVPIc:=proc(L,Ini,Fini,f,m,N)   local eq,var,a,i,Ini1, Li:

Ini1:=[op(Ini),seq(a[i],i=1..nops(Fini))]:

var:={seq(a[i],i=1..nops(Fini))}:

Li:=SIrecS(L,Ini1,f,m,N):

#Li[N+1]=Fini[1], Li[N]=Fini[2], Li[N-1]=Fini[3]:
eq:={seq(Li[N+2-i]=Fini[i], i=1..nops(Fini))}:

var:=solve(eq,var):

subs(var,Li):

end : 


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


#P1dClever(f,x,h,mu0,mu1) for cleverly solving Poisson 1-D equation
#u(x+h) = 2*u(x)-u(x-h)+u''(x) 

P1dClever:=proc(f,x,h,mu0,mu1)    local  eq,var,a,i,Ini1,Li   : 

[seq(dBVPIc([2,-1],[mu0],[mu1],2*f,i*x,1/h-1)[-1] , i=1..1/h) ]: 

end : 


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


#G1E(f,x,h,mu0,mu1) finds the global error in using mesh size h 

G1Eclever:=proc(f,x,h,mu0,mu1) local R1, A1clever:

R1:=P1ch(f,x,h,mu0,mu1):

A1clever:=P1dClever(f,x,h,mu0,mu1):

max(seq(abs(R1[i]-A1[i]),i=1..nops(R1))):

end:


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


#EmpiricalStability(f,x,mu0,mu1,ListH) empirically tests stability of G1Eclever

EmpiricalStability:= proc(f,x,mu0,mu1,ListH)   local h     : 

seq( G1Eclever(f,x,h,mu0,mu1)/h^2, h=op(ListH) ) : 

end: