%main
clear all;
clc;
fun = @sol_func;
x0 = [-0.1 -0.08 -0.05];
%x0 = [0,0,0];
P = fsolve(fun,x0);



function func = sol_func(P)

%Variable values
C1 = 1.07;                                              %F/m^2
C2 = 0.2;                                               %F/m^2
F = 96487;                                              %C/mol
R = 8.314;                                              %J/(K*mol)
T = 298;                                                %K
Ep = 80*8.854*10^-12;                                   %F/m
Ka1 = 1000*10^1.9;                                      %mol/m^3
Ka2 = 1000*10^-6.73;                                    %mol/m^3
Km = 0.001*10^-0.25;                                    %mol/m^3
Ntot = 8.3056*exp(-6);                                  %mol/m^2
pH = 6.5;
Ckci = .1;
                                                        %Equations

Hpb = 10^(-pH+3);                                       %(19)

    if(pH<=7)
        Kpb=Ckci;                                       %(20)
    else
        Kpb=Ckci-10^(-pH+3)+10^(-(14-pH)+3);            %(20)
    end

C0 = Hpb+Kpb;                                           %(18)

Kp = (Kpb*exp((-F/(R*T))*P(2)));                          %(4)
Hp = (Hpb*exp((-F/(R*T))*P(1)));                          %(5)


SiOH = (Ntot/(1+(Ka2/Hp)+(Hp/Ka1)+Km*Ka2*(Kp/Hp)));     %(11)

SiOHp2=(SiOH*Hp/Ka1);                                   %(7)
SiOm=(Ka2*SiOH/Hp);                                     %(8)
SiOK=(Km*(Ka2*Kp/Hp)*SiOH);                             %(9)

sig0 = F*(SiOHp2-SiOm-SiOK);                            %(15)
sigB = F*(SiOK);                                        %(16)
sigD = (-sign(P(3))*sqrt(2*R*T*Ep*C0*(exp(-F*...
    (P(3)/(R*T)))+exp(F*(P(3))/(R*T))-2)));             %(17)

% P(1) is psi0
% P(2) is psiB
% P(3) is psiD

func(1)= sig0 + sigB + sigD;
func(2)= sig0 - C1*(P(1) - P(2));
func(3)= sigD + C2*(P(2) - P(3));