clear all
close all

D = 1/6;
i = 10;  % mesh points
x = linspace(0,1,i);
dx = x(2)-x(1);  % mesh increment in meters
dt = 0.02;  % time step in seconds
n = 200; % number of time steps

f_n_minus_1_i = sin(2*pi*x);  % wave function starting condition t=0
% f_n_minus_1_i(end) = 0;
f_n_i = exp(-4*pi^2*D*dt)*sin(2*pi*x);  % wave function starting condition t=1*dt
% f_n_i(end) = 0;
f_n_i_right = [f_n_i(2:end) 0];  % wave function starting condition t=1*dt
f_n_i_left = [0 f_n_i(1:end-1)];  % wave function starting condition t=1*dt


% figure
% plot(f_n_minus_1_i)
% hold
% plot(f_n_i,'r')
% plot(f_n_i_right,'g')
% plot(f_n_i_left,'y')

for l = 2:(50/dt)

    f_n_plus_1_i = f_n_minus_1_i +D*2*dt/(dx^2)*(-2*f_n_i + f_n_i_right + f_n_i_left);
    f_n_plus_1_i(1) = 0;  % implementing BC
    f_n_plus_1_i(end) = 0;
    f_n_minus_1_i = f_n_i; %sets f old = f new
    f_n_i = f_n_plus_1_i;
    f_n_i_right = [f_n_plus_1_i(2:end) 0];
    f_n_i_left = [0 f_n_plus_1_i(1:end-1)];

    if (l ==0.06/dt)
        exact_sol_006 = exp(-4*pi^2*D*3*dt)*sin(2*pi*x);
        figure
        plot(f_n_i)
        hold
        plot(exact_sol_006,'r');
        LEGEND('computed','exact solution');
        title(['f as a function of x at t = ',num2str(0.06) ' seconds'])
        f_n_i_006 = f_n_i;

    end

    if (l ==0.1/dt)
        exact_sol_01 = exp(-4*pi^2*D*5*dt)*sin(2*pi*x);
        figure
        plot(f_n_i)
        hold
        plot(exact_sol_01,'r');
        LEGEND('computed','exact solution');
        title(['f as a function of x at t = ',num2str(0.1) ' seconds'])
        f_n_i_01 = f_n_i;

    end

    if (l ==0.9/dt)
        exact_sol_09 = exp(-4*pi^2*D*45*dt)*sin(2*pi*x);
        figure
        plot(f_n_i)
        hold
        plot(exact_sol_09,'r');
        LEGEND('computed','exact solution');
        title(['f as a function of x at t = ',num2str(0.9) ' seconds'])
        f_n_i_09 = f_n_i;

    end

    if (l ==50/dt)
        exact_sol_50 = exp(-4*pi^2*D*2500*dt)*sin(2*pi*x);
        figure
        plot(f_n_i)
        hold
        plot(exact_sol_50,'r');
        LEGEND('computed','exact solution');
        title(['f as a function of x at t = ',num2str(50) ' seconds'])
        f_n_i_50 = f_n_i;

    end

end