clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Use Galerkin method to approximately     %%%
%%% solve BVP -(exp(x)u')' = x, 0 < x < ell, %%%
%%%            u(0) = u(ell) = 0,            %%%
%%% on the n dimensional subspace            %%%
%%% Vn = span{phi_1,...,phi_n}, where        %%%
%%% phi_k(x) = sin(k*pi*x/ell).              %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = 500;  %%% dimension of subspace Vn
ell = 2;  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% compute the stiffness matrix in two ways %%%
%%% and compare how long each method took    %%%
%%% with tic and toc                         %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% First: the easy/naive way: %%%%%%%%%%%%%%%%%
%tic
%
%a = @(i,j,ell) i*j*pi^2*(ell^2+pi^2*i^2+pi^2*j^2)*(-1+(-1)^(i+j)*exp(ell)) / ( (ell^2 + pi^2*(i-j)^2)*(ell^2 + pi^2*(i+j)^2) );
%
%for i = 1:n
%    for j = 1:n
%        K(i,j) = a(i,j,ell);
%    end
%end
%toc
%
%%% The cleverer way:    %%%%%%%%%%%%%%%%%%%%%%%
tic
aa = @(i,j,ell) (i*j*pi^2).*((ell^2+pi^2*(i.^2*ones(1,n)+ones(n,1)*j.^2)).*(-1+(-1).^(i*ones(1,n)+ones(n,1)*j)*exp(ell))) ./ ( (ell^2 + pi^2*(i*ones(1,n)-ones(n,1)*j).^2).*(ell^2 + pi^2*(i*ones(1,n)+ones(n,1)*j).^2) );

K = aa((1:n)',(1:n),ell);
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% fcoeffs(i) = (f,phi_i) [L^2 inner product] %
%%% where f(x) = x.                            %
%%% Work out on paper, or using Maple,         %
%%% Mathematica, or symbolic toolobox.         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fcoeffs = @(i) (-1).^(i+1)*ell^2./i/pi;
f = fcoeffs(1:n)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Solve  K*c = fcoeffs  %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c =  K\f;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% calculate Galerkin approx               %%%%
%%% un = c(1)phi_1(x) + ... + c(1)phi_1(x)  %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x = linspace(0,ell,1000)';
un = 0*x;
for k = 1:n
    un = un + c(k)*sin(k*pi*x/ell);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% uactual = actual solution evaluated at gridpts x. %
%%% Work out on paper, or using Maple,                %
%%% Mathematica, or symbolic toolobox.                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
uactual = (1/2)*exp(-x).*(-(2*exp(-ell)-2+2*exp(-ell)*ell+exp(-ell)*ell^2)/(exp(-ell)-1)+2+2*x+x.^2)+(1/2)*exp(-ell)*ell*(2+ell)/(exp(-ell)-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1), clf, box on, hold on
plot(x,uactual,'r:','LineWidth',2)
plot(x,un','k')
legend('solution','approx')
xlabel('x','fontsize',15)
title(sprintf('n = %d',n),'fontsize',15)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2), clf, box on, hold on
plot(x,abs(uactual-un))
xlabel('x','fontsize',15)
title(sprintf('error |u-un| for n = %d',n),'fontsize',15)