% input: mu - relation between space step h and time step k given by 
%             mu = k*h^-2
%        h  - space step
%        T  - final time when the solution is plotted
% output: v - finite difference solution at time T

function v = cn(mu,h,T)

b =1;


k = mu*h^2; % set the number oftime step
len = 10; % length of the interval of interest
m = len/h; % number of space grid points
n = floor(T/k); % number of time grid points

% initial conditions: square function
vint = zeros(m-1,1);
for i = floor(m/2-m/len):floor(m/2+m/len)
    vint(i) = 1;
end

% boundary conditions
data1 = zeros(n,1);
data2 = zeros(n,1);
% bounday condition of the exact solution
for i = 1:n
data1(i) = .5*(erf((1+len/2)/sqrt(4*i*k)) + ...
    erf((1-len/2)/sqrt(4*i*k)));
data2(i) = .5*(erf((1+len/2)/sqrt(4*i*k)) + ...
    erf((1-len/2)/sqrt(4*i*k)));
end

v = vint;    

% initialization for Thomas Algorithm
p = zeros(m,1);
q = zeros(m,1);
aa = -b*mu/2;
bb = 1+b*mu;
cc = aa;

for i = 1:n  
    % p and q on left boundary
    p(1) = 0;
    q(1) = data1(i);
    dd = v(1) + .5*b*mu*(v(2) - 2*v(1) + data1(i));
    denom = aa*p(1) + bb;
    p(2) = -cc/denom;
    q(2) = (dd - q(1)*aa)/denom;
    
    % p and q on right boundary
    dd = v(m-1) + .5*b*mu*(data2(i) - 2*v(m-1) + v(m-2));
    denom = aa*p(m-1) + bb;
    p(m) = -cc/denom;
    q(m) = (dd - q(m-1)*aa)/denom;
    
    % p and q on inner grid points
    for j = 2:m-2
        dd = v(j) + .5*b*mu*(v(j+1) - 2*v(j) + v(j-1));
        denom = aa*p(j) + bb;
        p(j+1) = -cc/denom;
        q(j+1) = (dd - q(j)*aa)/denom;
    end
    
    % function value at right boundary
    v(m-1) = p(m)*data2(i) + q(m);
    
    % function value atinner grid points
    for k=m-2:-1:1
        v(k) = p(k+1)*v(k+1) + q(k+1);
    end
end

% excat solution
vexact = zeros(m-1,1);
for i = 1:m-1
    vexact(i) = .5*(erf((1-(i*h-len/2))/sqrt(4*T)) + ...
        erf((1+(i*h-len/2))/sqrt(4*T)));
end

% plotting
x = (-len/2:h:len/2)';
v = [data1(n) ; v ; data2(n)];
vint = [data1(n) ; vint ; data2(n)];
vexact = [data1(n) ; vexact ; data2(n)];
figure(2)
subplot(2,1,1)
hold on
title('Initial condition and fin. diff. solution at time T for the Crank-Nicolson scheme')
plot(x,v,'r')
plot(x,vint,'b')
hold off
legend('fin. diff. soln', 'initial condition')
subplot(2,1,2)
plot(x,abs(vexact-v),'k')
legend('error')
hold off