% parameter values
lambda = 0.1;
alpha  = 1;
sigma  = 0.2;

% nodal points
SI = [1/3;2/3];
SB = [0;1];

% basis matrices
PhiI0   = [ones(size(SI,1),1) SI SI.^2 SI.^3];
PhiB0   = [ones(size(SB,1),1) SB SB.^2 SB.^3];
PhiI1   = [zeros(size(SI,1),1) ones(size(SI,1),1) 2*SI 3*SI.^2];
PhiI11  = [zeros(size(SI,1),2) 2*ones(size(SI,1),1) 6*SI];

% compute approximation coefficients
f = ones(size(SI,1),1);
g = SB;
L = diag(SI)*PhiI0 + diag(lambda*(alpha-SI))*PhiI1+sigma*PhiI11;
c = [L;PhiB0]\[f;g];

% evaluate function and residual
s      = linspace(0,1,101)';
PhiI0  = [ones(size(s,1),1) s s.^2 s.^3];
PhiI1  = [zeros(size(s,1),1) ones(size(s,1),1) 2*s 3*s.^2];
PhiI11 = [zeros(size(s,1),2) 2*ones(size(s,1),1) 6*s];
L      = diag(s)*PhiI0 + diag(lambda*(alpha-s))*PhiI1+sigma*PhiI11;
r      = L*c-1;

% Plot solution approximation and residual function
figure(1)
plot(s,PhiI0*c,'b');
title('Appproximation Solution')
xlabel('S')
ylabel('V(S)')

figure(2)
plot(s,r,'b');
title('Residual Function')
xlabel('S')
ylabel('r(S)')


% Compute higer order approximation
fspace=fundefn('spli',21,0,1);
s=funnode(fspace);
SI=s(2:end-1);
SB=[0;1];
PhiI0=funbas(fspace,SI);
PhiB0=funbas(fspace,SB);
PhiI1=funbas(fspace,SI,1);
PhiI11=funbas(fspace,SI,2);

f=ones(size(SI,1),1);
g=SB;
L=diag(SI)*PhiI0 + diag(lambda*(alpha-SI))*PhiI1+sigma*PhiI11;
c=[L;PhiB0]\[f;g];

s=linspace(0,1,101)';
PhiI0=funbas(fspace,s);
PhiI1=funbas(fspace,s,1);
PhiI11=funbas(fspace,s,2);
L=diag(s)*PhiI0 + diag(lambda*(alpha-s))*PhiI1+sigma*PhiI11;
r=L*c-1;

figure(1)
hold on
plot(s,PhiI0*c,'r');
hold off
title('Appproximation Solution')
xlabel('S')
ylabel('V(S)')

figure(2)
hold on
plot(s,r,'r');
hold off
title('Residual Function')
xlabel('S')
ylabel('r(S)')