#hw23.txt, 19 April 2014, Andrew Lohr 

Pw:=proc(Lf,x,Lx) local i,command:
command:= x<Lx[1],0:
for i from 2 to nops(Lx) do
command := command,x>=Lx[i-1] and x<Lx[i],Lf[i-1]:
od:
piecewise(command,x>=Lx[nops(Lx)],0):
end: 

Compare :=proc(p,q,v,x,x1) local u,f,N,ret:
u := v*x*(1-x):
f:=-diff(p*diff(u,x),x) + q*u:

ret:=[]:
for N in {10,20,30,40} do
ret:=[op(ret),RR(p,q,f,x,[seq(i/N,i=0..N)])[round(x1*N+1)]]: 
od:
evalf(subs(x=x1,[op(ret),u])):
end:

#running it with 
#    Compare(x^2,1+x^3,sin(x),x,.513);
#    Compare(x^2,1+x^3,x^7+1,x,.513);
#    Compare(x^2+3,sin(x),exp(x),x,.542); 

#we get
#   [0.1251017497, 0.1234668473, 0.1231112603, 0.1227372190, 0.1226153759]
#   [0.2562016786, 0.2535615542, 0.2529879707, 0.2522113970, 0.2521669733]
#   [0.4355599517, 0.4264488333, 0.4274669885, 0.4266447327, 0.4268274813]


BSpline:=proc(d,Lx,i) local a,Lu,i1,j1,var,eq:
if nops(Lx)<>d+2 then
  RETURN(FAIL):
fi:
Lu:=[seq(add(a[i1,j1]*x^j1,j1=0..d),i1=0..d)]:
var:={seq(seq(a[i1,j1],j1=0..d),i1=0..d)}:
eq:={subs(x=Lx[1],Lu[1]),
seq(subs(x=Lx[1],diff(Lu[1],x$i1)),i1=1..d-1)}:
eq:= eq union {subs(x=Lx[d+2],Lu[d+1]),
seq(subs(x=Lx[d+2],diff(Lu[d+1],x$i1)),i1=1..d-1)}:
#must be 1 at proscribed position
eq:= eq union {subs(x=Lx[i+1],Lu[i] = 1)}:
#smoothness conditions
eq:= eq union {seq(subs(x=Lx[i1],Lu[i1] = Lu[i1-1]),i1 = 2..d+1)}:
eq := eq union
	{seq( 
	seq( subs(x=Lx[i1],diff(Lu[i1-1],x$j1)) - subs(x=Lx[i1],diff(Lu[i1],x$j1)),j1=1..d-1)
	,i1= 2..d+1)}:
var:=solve(eq,var):
if(var = NULL) then
return(FAIL):
fi:
Pw(subs(var,Lu),x,Lx):
end:

TestF3:=proc(c,Lx,x) local i:
     expand(add(c[i]*BSpline(nops(Lx)-2,Lx,i+1),i=1..nops(Lx)-2)), 
      { seq( c[i], i=1..nops(Lx)-2)}:
     end:

 RR3:=proc(p,q,f,x,Lx) local c, Lu, F, eq,var:
    Lu:=TestF3(c,Lx,x):
    var:=Lu[2]:
    Lu:=Lu[1]: 
   F:=EvalFunct(p,q,f,x,Lx,Lu):
    eq:=evalf({seq(diff(F,var[i]),i=1..nops(var))}):   
    var:=solve(eq,var):
   if var=NULL then
       FAIL:
     else
       subs(var,Lu):
     fi:  
   end: