% 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 = ftcs(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

vold = vint;    
vnew = zeros(m-1,1);

% loop over time and space
for i=1:n
    for j = 2:m-1
        
        % function value at left boundary
        if j == 2
            vnew(j) = (1-2*b*mu)*vold(j) +b*mu*(vold(j+1)-data1(i));
        % function value at right boundary
        elseif j == m-1
            vnew(j) = (1-2*b*mu)*vold(j) +b*mu*(vold(j-1)-data2(i));
        % function value at inner grid points
        else
            vnew(j) = (1-2*b*mu)*vold(j) +b*mu*(vold(j+1)+vold(j-1));
        end
        
    end
    vold = vnew;
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) ; vold ; 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 forward-time central-space 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