/* interval */
a:0$ b:1$ [a,b];

/* subdivisions */
n:5;
dx:(b-a)/n;

/* 3 points give the footprint of each tent */
bb(i):=a+i*dx$ /* center point */
aa(i):=bb(i)-dx$
cc(i):=bb(i)+dx$

/* shape functions */
bf(i,x):=(x-aa(i))*charfun(aa(i)<x and x<bb(i))
    -(x-cc(i))*charfun(bb(i)<x and x<cc(i));

/* [ characteristic function:
    charfun(true)=1, charfun(false)=0 ] */

create_list(bf(i,x),i,1,n-1)$
plot2d(%,[x,a,b])$

/* w -> <w,bf(j)> */
gi(j,w):=integrate(w*(x-aa(j)),x,aa(j),bb(j))
    -integrate(w*(x-cc(j)),x,bb(j),cc(j));

/* w -> <w,bf(j)'> */
gi1(j,w):=integrate(w,x,aa(j),bb(j))
    -integrate(w,x,bb(j),cc(j));

/* arbitrary v in span{shape functions} */
v(x):=sum(c[i]*bf(i,x),i,1,n-1);

/* residual for v''+v=1
0=<R(v),bf(j)>=<v'',bf(j)>=-<v',bf(j)'>+<v+1,bf(j)> */
eqs:create_list(-gi1(j,diff(v(x),x))+gi(j,v(x)+1),j,1,n-1);

/* approximate solution */
vars:create_list(c[i],i,1,n-1);
solve(eqs,vars);
as:subst(%,v(x));

/* exact solution */
es:cos(x)+(1-cos(1))/sin(1)*sin(x)-1;

plot2d([as,es],[x,a,b])$

/* coefficient matrix -<bf(i)',bf(j)'>+<bf(i),bf(i)> */
/* augmentation column: -<1,bf(j)> */
genmatrix(lambda([i,j],-gi1(j,diff(bf(i,x),x))+gi(j,bf(i,x))),n-1,n-1)$
A:addcol(%,transpose(genmatrix(lambda([i,j],-gi(j,1)),1,n-1)));

/* reduced row echelon form */
rref(a):=block([p,q,k],[p,q]:matrix_size(a),a:echelon(a),
    k:min(p,q),
    for i thru min(p,q) do (if a[i,i]=0 then (k:i-1,return())),
    for i:k thru 2 step -1 do (for j from i-1 thru 1 step -1 do a:rowop(a,j,i,a[j,i])),
    a)$
rref(A);