%Starting the script by using the hints from the hint sheet given in the
%exercise hours

close all; clear; clc

Nz = 150;
dx = 1 /(Nz-1);
[x,w,A,B] = FDGrid(Nz,dx);
PAR1 = 100; PAR2 = 1000;
xA0 = ones(Nz,1);
xA = fsolve(@(xA) CPB(Nz, xA, A, B, PAR1, PAR2), xA0);

xB = 1-xA;


disp(xA)

FigProps = struct('Color',[1,1,1],'OuterPosition',[120,120,1280,690]);
FontProps = struct('FontName', 'Calibri', 'FontSize', 18, 'FontWeight','bold');
fig = figure(1);
set(fig, FigProps);
set(axes, FontProps);
set(plot(x,xA), 'Color', 'r', 'LineWidth', 2, 'DisplayName','x_A');hold on;
set(plot(x,xB), 'Color', 'r', 'LineWidth', 2,'Marker','*', 'DisplayName','x_B');hold on;
set(xlabel('z/L'),FontProps);set(ylabel('Mole fraction'),FontProps);
set(title('Mole fractions xA and xB from CPB Model, PAR1 = 100; PAR2 = 1000'),FontProps);
xlim([0 1]);ylim([0 1]); box on; grid on;
legend




function [f] = CPB(Nz, xA, A, B, PAR1, PAR2)

    d1 = A*xA; % This is the first order derivative
    d2 = B*xA; % This is the second order derivative
    
    f(1) = xA(1)-1;
    
    for i = 2:(Nz-1)
       f(i) = d2(i) - PAR1*d1(i) - PAR2*log(1/(1-0.5*xA(i)));
    end
    
    f(Nz) = d1(Nz);

end