%This script numerically solves the diffusion equation and plots the
%results.

clear

m = input('Please enter the number of lines (>= 3): ');
m = floor(m);
while(m < 3)
    m = input('Please enter the number of lines (>= 3): ');
end

tf = input('Please enter the final time (> 0): ');
while(tf <= 0)
    tf = input('Please enter the final time (> 0): ');
end

tNum = input('Please enter the number of divisions of the time interval (>= 1): ');
tNum = floor(tNum);
while(tNum < 1)
    tNum = input('Please enter the number of divisions of the time interval (>= 1): ');
end

dx = 1 / (m + 1);
initC = zeros(m, 1);
for i = 1:m
    initC(i) = initialCondition(i * dx);
end

[t, y] = ode45('discreteVector', 0:(tf / tNum):tf, initC);

plotX = zeros((m + 2), (tNum + 1));
plotY = zeros((m + 2), (tNum + 1));
plotZ = zeros((m + 2), (tNum + 1));
dx = 1 / (m + 1);
dt = 1 / tNum;

for i = 1:(tNum + 1)
    plotX(1, i) = t(i);
    plotY(1, i) = 0;
    plotZ(1, i) = leftBoundary(t(i));
    plotX((m + 2), i) = t(i);
    plotY((m + 2), i) = 1;
    plotZ((m + 2), i) = rightBoundary(t(i));
    for j = 2:(m + 1)
        plotX(j, i) = t(i);
        plotY(j, i) = (j - 1) * dx;
        plotZ(j, i) = y(i, (j - 1));
    end
end

plot3(plotX, plotY, plotZ);