clear;
close all;
clc;

% Task 2

% Creating grid
N = 200;
[z_hat, w, A, B] = FiniteDifferenceGrid (N);

% Initital guess
xA0 = 0.5*ones(N,1);

% Solving for constant alpha, while varying beta:
xA1 = fsolve(@(xA1) sysEquations(xA1, A, B, 100, 10, N), xA0);
xA2 = fsolve(@(xA2) sysEquations(xA2, A, B, 100, 100, N), xA0);
xA3 = fsolve(@(xA3) sysEquations(xA3, A, B, 100, 500, N), xA0);
xA4 = fsolve(@(xA4) sysEquations(xA4, A, B, 100, 1000, N), xA0);

xB1 = 1-xA1;
xB2 = 1-xA2;
xB3 = 1-xA3;
xB4 = 1-xA4;


%% Plotting the result

% Properties of the plot
figProps = struct('Color', [1 1 1], 'OuterPosition', [170, 170, 1000, 700]);
fontProps = struct('FontName', 'Calibri', 'FontSize', 18, 'FontWeight', 'bold');
% Initializing the plot
fig = figure(1);
ax = axes;
% Plotting the results
plot(z_hat, xA1, 'Color', 'r', 'LineWidth', 2, 'Marker', 'o', 'MarkerEdgeColor', 'r','Markersize', 4, 'MarkerFaceColor', 'none', 'DisplayName', 'x_A, \beta = 10');
hold on;
plot(z_hat, xB1, 'Color', 'r', 'LineWidth', 2, 'Marker', 'none', 'DisplayName', 'x_B, \beta = 10');
plot(z_hat, xA2, 'Color', "#EDB120", 'LineWidth', 2, 'Marker', 'o', 'MarkerEdgeColor', "#EDB120",'Markersize', 4, 'MarkerFaceColor', 'none', 'DisplayName', 'x_A, \beta = 100');
plot(z_hat, xB2, 'Color', "#EDB120", 'LineWidth', 2, 'Marker', 'none', 'DisplayName', 'x_B, \beta = 100');
plot(z_hat, xA3, 'Color', "#77AC30", 'LineWidth', 2, 'Marker', 'o', 'MarkerEdgeColor', "#77AC30",'Markersize', 4, 'MarkerFaceColor', 'none', 'DisplayName', 'x_A, \beta = 500');
plot(z_hat, xB3, 'Color', "#77AC30", 'LineWidth', 2, 'Marker', 'none', 'DisplayName', 'x_B, \beta = 500');
plot(z_hat, xA4, 'Color', "#4DBEEE", 'LineWidth', 2, 'Marker', 'o', 'MarkerEdgeColor', "#4DBEEE",'Markersize', 4, 'MarkerFaceColor', 'none', 'DisplayName', 'x_A, \beta = 1000');
plot(z_hat, xB4, 'Color', "#4DBEEE", 'LineWidth', 2, 'Marker', 'none', 'DisplayName', 'x_B, \beta = 1000');
% Tidying up
set(fig, figProps);
set(ax, 'FontSize', fontProps.FontSize);
set(xlabel("z' = z/L", 'interpreter', 'tex'), fontProps);
set(ylabel('Molar fractions, x_i'), fontProps);
set(title('Molar fractions x_A and x_B, while varying \beta, and keeping \alpha = 100','interpreter', 'tex'), fontProps);
xlim([0 1]);
box("on");
grid("on");
set(legend(), 'FontName', fontProps.FontName, 'interpreter', 'tex', 'Location', 'northeastoutside')
%Saving the figure
saveas(fig, 'Ex6_3_plot_1', 'jpg')


% Solving for constant beta, while varying alpha:
xA5 = fsolve(@(xA5) sysEquations(xA5, A, B, 10, 1000, N), xA0);
xA6 = fsolve(@(xA6) sysEquations(xA6, A, B, 50, 1000, N), xA0);
xA7 = fsolve(@(xA7) sysEquations(xA7, A, B, 100, 1000, N), xA0);
xA8 = fsolve(@(xA8) sysEquations(xA8, A, B, 200, 1000, N), xA0);

xB5 = 1-xA5;
xB6 = 1-xA6;
xB7 = 1-xA7;
xB8 = 1-xA8;


%% Plotting the result

% Initializing the plot
fig = figure(2);
ax = axes;
% Plotting the results
plot(z_hat, xA5, 'Color', 'r', 'LineWidth', 2, 'Marker', 'o', 'MarkerEdgeColor', 'r','Markersize', 4, 'MarkerFaceColor', 'none', 'DisplayName', 'x_A, \alpha = 10');
hold on;
plot(z_hat, xB5, 'Color', 'r', 'LineWidth', 2, 'Marker', 'none', 'DisplayName', 'x_B, \alpha = 10');
plot(z_hat, xA6, 'Color', "#EDB120", 'LineWidth', 2, 'Marker', 'o', 'MarkerEdgeColor', "#EDB120",'Markersize', 4, 'MarkerFaceColor', 'none', 'DisplayName', 'x_A, \alpha = 100');
plot(z_hat, xB6, 'Color', "#EDB120", 'LineWidth', 2, 'Marker', 'none', 'DisplayName', 'x_B, \alpha = 100');
plot(z_hat, xA7, 'Color', "#77AC30", 'LineWidth', 2, 'Marker', 'o', 'MarkerEdgeColor', "#77AC30",'Markersize', 4, 'MarkerFaceColor', 'none', 'DisplayName', 'x_A, \alpha = 500');
plot(z_hat, xB7, 'Color', "#77AC30", 'LineWidth', 2, 'Marker', 'none', 'DisplayName', 'x_B, \alpha = 500');
plot(z_hat, xA8, 'Color', "#4DBEEE", 'LineWidth', 2, 'Marker', 'o', 'MarkerEdgeColor', "#4DBEEE",'Markersize', 4, 'MarkerFaceColor', 'none', 'DisplayName', 'x_A, \alpha = 1000');
plot(z_hat, xB8, 'Color', "#4DBEEE", 'LineWidth', 2, 'Marker', 'none', 'DisplayName', 'x_B, \alpha = 1000');
% Tidying up
set(fig, figProps);
set(ax, 'FontSize', fontProps.FontSize);
set(xlabel("z' = z/L", 'interpreter', 'tex'), fontProps);
set(ylabel('Molar fractions, x_i'), fontProps);
set(title('Molar fractions x_A and x_B, while varying \alpha, and keeping \beta = 1000','interpreter', 'tex'), fontProps);
xlim([0 1]);
box("on");
grid("on");
set(legend(), 'FontName', fontProps.FontName, 'interpreter', 'tex', 'Location', 'northeastoutside')
%Saving the figure
saveas(fig, 'Ex6_3_plot_2', 'jpg')


%% Functions used in the problem

% Creating a function for the f(xA) defined on my submission
function toReturn = f(x)
   toReturn = log(1/(1-0.5*x)); 
end

% Creating the equation system to be solved
function wantToBeZero = sysEquations(xA, dxA_matrix, ...
    d2xA2_matrix, param1, param2, numberOfNodes)
    
    %Creating empty array to optimize performance
    wantToBeZero = zeros(numberOfNodes, 1);

    % Computing the differentials values beforehand using the grid
    d2xA_dz2 = d2xA2_matrix*xA;
    dxA_dz = dxA_matrix*xA;

    %% Boundary conditions
    % At z=0
    wantToBeZero(1) = xA(1)-1;

    %% General solution for the middle nodes:
    for i = 2:(numberOfNodes-1)
        wantToBeZero(i) = d2xA_dz2(i) - param1*dxA_dz(i) - param2*f(xA(i));
    end
    % At z=1
    wantToBeZero(numberOfNodes) = dxA_dz(numberOfNodes);
end 