% Program for solving problem of propagation of a time-harmonic plane wave
%   along the x-axis
%
% -u''-k^2 u = 1, x\in[0,1], u(0)=0, u'(1)-iku(1)=0

k=input('Enter value of k:- '); %enter wavenumber
N=input('Enter value of N:- '); %enter number of elements


h=1/N;   % set mesh width = 1/N
A=sparse((2/h)*eye(N)-(1/h)*diag(ones(N-1,1),1)-(1/h)*diag(ones(N-1,1),-1));
A(N,N)=1/h; % compute A := \int_0^1 \chi_j'(x) \chi_m'(x) dx
B=sparse((2*h/3)*eye(N)+(h/6)*diag(ones(N-1,1),1)+(h/6)*diag(ones(N-1,1),-1));
B(N,N)=h/3; % compute B := \int_0^1 \chi_j(x) \chi_m(x) dx
C=sparse(zeros(N));C(N,N)=1;  %compute C:=u(1)*\chi_m(1), from boundary data

LHS=(A-k^2*B-i*k*C); %put them all together to get the matrix on the left hand side.
RHS=h*ones(N,1);RHS(N)=h/2; %evaluate RHS := \int_0^1 \chi_j(x) dx

ucoeff=LHS\RHS; % Solve LHS*ucoeff=RHS, to get coefficients for solution.

% Compute true solution at 10*k points

x=linspace(0,1,10*k);
utrue=(exp(i*k*x)-1-i*exp(i*k)*sin(k*x))/(k^2);

% The rest is for computing U(x) given the coefficients U_1, U_2, etc

grid=0:h:1;
uapprox=zeros(size(x));
for j=1:length(x);
        r=max(1,sum(x(j)>grid));
    if r==1;
        uapprox(j)=ucoeff(1)*x(j)/h;
    else
        uapprox(j)=(ucoeff(r)*(x(j)-grid(r))/h)+...
            (ucoeff(r-1)*(grid(r+1)-x(j))/h);
    end
end
    
% Plot some results

subplot(221)
plot(x,real(utrue),':',x,real(uapprox));legend('u','U','location','SouthEast');...
    title('real part of true and approx')
subplot(222)
plot(x,imag(utrue),':',x,imag(uapprox));legend('u','U','location','SouthEast');...
    title('imaginary part of true and approx')
subplot(212)
plot(x,abs(utrue-uapprox)/max(abs(utrue)));title('Relative Error |u-U|/max|u|')