The PDE to be solved is:
 Uxx(x) + 2*Ux(x) + U(x)^2 = f(x)
 f(x) = x^4 + 2x^2 + 4x + 3
 discrete u(xmin=0)=1  u(xmax=1)=2
 analytic boundary ub(x) = x^2+1
 
xg[0]=0.0, uA[0]=1.0
xg[1]=0.25, uA[1]=1.0625
xg[2]=0.5, uA[2]=1.25
xg[3]=0.75, uA[3]=1.5625
xg[4]=1.0, uA[4]=2.0
 
i=1, ia=0, x=0.25
A[0][0]=-33.33333333333333
A[0][1]=20.0
A[0][2]=1.333333333333333
A[0][3]=1.0
A[0][4]=0.0
A[0][5]=0.0
Y[0]=-7.204427083333333
 
i=2, ia=1, x=0.5
A[1][0]=16.0
A[1][1]=-40.0
A[1][2]=26.666666666666664
A[1][3]=0.0
A[1][4]=1.0
A[1][5]=0.0
Y[1]=10.229166666666668
 
i=3, ia=2, x=0.75
A[2][0]=9.333333333333332
A[2][1]=-4.0
A[2][2]=-20.0
A[2][3]=0.0
A[2][4]=0.0
A[2][5]=1.0
Y[2]=-23.891927083333332
 
simeq_newton5 itr 1, prev=20.583333333333332, residual=0.34766823255646795
simeq_newton5 itr 2, prev=0.34766823255646795, residual=0.001161531636257429
simeq_newton5 itr 3, prev=0.001161531636257429, residual=7.264450374577791E-9
check solution against PDE
i=1, x=0.25, u=1.0625
 ux=0.4999999992766889, uxx= 2.0000000043552073
ubx=0.5, ubxx=2.0
uAx=0.5000000000000006, uAxx=1.9999999999999996
eval=4.128906252273183, f(x)=4.12890625
err=2.2731834192768474E-9
ederiv=4.128906250000001, f(x)=4.12890625
eerr=8.881784197001252E-16
aderiv=4.12890625, f(x)=4.12890625
aerr=0.0
 
i=2, x=0.5, u=1.25
 ux=1.0000000003556657, uxx= 2.0000000037627657
ubx=1.0, ubxx=2.0
uAx=0.9999999999999999, uAxx=1.9999999999999942
eval=5.562500003619391, f(x)=5.5625
err=3.6193910091242287E-9
ederiv=5.562499999999994, f(x)=5.5625
eerr=-6.217248937900877E-15
aderiv=5.5625, f(x)=5.5625
aerr=0.0
 
i=3, x=0.75, u=1.5625
 ux=1.500000000901144, uxx= 2.0000000000872156
ubx=1.5, ubxx=2.0
uAx=1.5, uAxx=2.0
eval=7.441406251371886, f(x)=7.44140625
err=1.3718857161393316E-9
ederiv=7.44140625, f(x)=7.44140625
eerr=0.0
aderiv=7.44140625, f(x)=7.44140625
aerr=0.0
 
maxerr=3.6193910091242287E-9, rmserr=2.5916162622805875E-9, avgerr=2.4214867148468025E-9