# 	Here, solve -u'' + p(x)u' + q(x)u = f, u(0) = a, u'(1) = b, 
# 	by first remapping the problem to be homogeneous, and using
# 	the sinc basis plus two splines in a collocation procedure.
#
# 
Digits := 10;

read `sinc8.mas`;		# all subroutines

w_1 := unapply(poly(0,1,0,0, 0,1), x);
w_2 := unapply(poly(0,0,1,0, 0,1), x);

# 	The shift:
a := 2; 
b := 1;		# choose various cases for testing
Q := poly(1,0,0,1,0,1);
P := poly(0,0,0,1,0,1);
PHI := a*Q + (b-a)*P;

phi := x -> ln(x/(1-x));		# from [0,1] to R
z = phi(u):
solve(", u):
psi := unapply(", z);			# from R to [0,1]

# g1 := unapply(1/diff(phi(x),x),x);
g1 := unapply(x**(9/10)*(1-x)**(9/10), x);
Sphi := (i,h,x) -> sinc((phi(x)-i*h)/h);

define(L, forall([g], L(g) = -Diff(g,x,x) + p(x)*Diff(g,x) + q(x)*g));
vv := proc(g);		# value(value(arg))
   value(value(g));
end;

#
p   := x -> -1/(6*x);		# the functions in the D.E.
q   := x -> 1/x**2;
f   := x -> -1/12*(-57*x^3+66*x^4+21*x^(3/2)*a-10*x^(3/2)*(1-x)^(1/2)*b-12*x^(1/2)*a)/x^(5/2)/(1-x)^(1/2);
exa := x -> (a+x^(3/2))*(1-x)^(3/2)+b*x;	# exact answer, for checking
check_should_be_zero := factor(vv(L(exa(x)) - f(x)));

f(x) - vv(L(PHI)):
f := unapply(", x);

# 	system sizes:
M := 7;
N := 7;
h := evalf(Pi/sqrt(2*M));
m := M+N+1;

# 	Define the "coefficient functions" of the actual coefficients:
# 
k := 'k';
j := 'j';	# clear for the following definitions:
coe[spline1] := vv(L(w_1(x)))*Sphi(j,h,x);
coe[sincs] := vv(L(sinc((phi(x)-k*h)/h)*g1(x)))*Sphi(j,h,x);
coe[spline2] := vv(L(w_2(x)))*Sphi(j,h,x);

# 	now, get final lhs, rhs and solve system to get u
#	(this is the actual "integration")
with(linalg):
llhs := matrix(m+2, m+2, 0):
rrhs := matrix(m+2, 1, 0):
for j from -M-1 to N+1 do;	# inner products with Sphi(j,h,x)
				# (also the row index)
   indj := j+M+2;
   xj := evalf(psi(j*h));
   # 	the 1st spline:
   k := -M-1;
   indk := k+M+2;
   tmp := evalf(subs(x=xj, coe[spline1]));
   llhs[indj,indk] := tmp;
   #	the sinc series:
   for k from -M to N do;	# the coefficients (column index)
      indk := k+M+2;
      tmp := evalf(subs(x=xj, coe[sincs]));
      llhs[indj,indk] := tmp;
   od; 	# k
   k := N+1;
   indk := k+M+2;
   tmp := evalf(subs(x=xj, coe[spline2]));
   llhs[indj,indk] := tmp;
   # 	The right side:
   rrhs[indj,1] := evalf(f(xj)*Sphi(j,h,xj));
od;

#  	solve system:
mat := blockmatrix(1,2, llhs, rrhs):
gausselim(mat):
u_vec := backsub(");

#	and form the computed answer from components:
k := 'k';
sum('u_vec[k+M+2]*sinc((phi(x)-k*h)/h)*g1(x)', 'k = -M..N'):
u := unapply(evalf(" + PHI + u_vec[1]*w_1(x) + u_vec[m+2]*w_2(x)), x);
up := unapply(diff(u(x),x),x):

exa_p := unapply(diff(exa(x),x),x):


# 	approximate the original function:
p1 := plot(u(x),   x=0..1, style=point):
p2 := plot(exa(x), x=0..1, style=line):
plots[display]({p1,p2}, title = cat(`F,   N = `, N, `    Digits = `, Digits));

p1 := plot(up(x),   x=0..1, style=point):
p2 := plot(exa_p(x), x=0..1, style=line):
plots[display]({p1,p2}, title = cat(`F,   N = `, N, `    Digits = `, Digits));