function strongVariationalFormulation(n)

% Find u in H^2(0,1): u(0)=0, -u'(1)=1 and
% int_0^1 -u''(x) v(x) dx = \int_0^1 1 v(x) dx for all v in L^2(0,1).

u = @(x) (-0.5*x.^2);

h = 1/n;

phi1 = @(x) (1-3*x.^2+2*x.^3);
phi2 = @(x) (3*x.^2-2*x.^3);
ksi1 = @(x) (x-2*x.^2+x.^3);
ksi2 = @(x) (-x.^2+x.^3);
psi0 = @(x) (x.^0);
psi1 = @(x) (2*x-1);

Ae = -(1/h) * [0 0 -1 1;2 -2 1 1];
be = h*[1;0];

Atilde = zeros(2*n,2*(n+1));
btilde = zeros(2*n,1);
for i=1:n
    idxu = [[i,i+1],n+1+[i,i+1]];
    idxv = [i,n+i];
    k = 1;
    Atilde(idxv,idxu) = Atilde(idxv,idxu) + k*Ae;
    btilde(idxv) = btilde(idxv) + be;
end

A = Atilde(:,2:(end-1));
b = btilde - (-1)*h*Atilde(:,end);

uuprime = [0;A\b;-1*h];

uh = uuprime(1:(end/2));
uprimeh = uuprime(((end/2)+1):end);

xx = 0:0.01:1;
hh = plot(xx,u(xx),'r',0:h:1,uh','x');
% set(hh,'LineWidth',5)
hold on
for i=1:n
    x = h*(i-1)+h*xx;
    hh = plot(x,uh(i)*phi1(xx)+uh(i+1)*phi2(xx)+uprimeh(i)*ksi1(xx)+uprimeh(i+1)*ksi2(xx),'b--');
    set(hh,'LineWidth',5)
end