Demo 10.4.2

u All the coefficients and response need to be functions:

exact = x -> exp(sin(x));
p = x -> -cos(x);
q = sin;
r = x -> 0; 

x,u = FNC.bvplin(p,q,r,[0,3*pi/2],1,exp(-1),30);

plot(exact,0,3π/2,layout=(2,1),label="exact");
scatter!(x,u,m=:o,subplot=1,label="numerical",
    yaxis=("solution"),title="Solution of a linear BVP");

plot!(x,exact.(x)-u,subplot=2,xaxis=L"x",yaxis=("error"))