function [b,B,x2,pH]=Newtonprotolys(x1)
%Newtons metod för en protolysjämvikt i vattenlösning. Givet
%gissningsvärden på koncentrationerna av HA, A, OH och H i (4x1)-matrisen x1 räknar programmet
%ut residualen, Jacobimatrisen A, en bättre gissning x2 och pH.

%Följande ekvationer ska uppfyllas:
% HA <-> A + H      Ka = [A]*[H]/[HA]       Syrans protolysjämvikt
% nHA0 + nA0 = V*([HA] + [A])               Massbalans 
% H2O <-> OH+H    Kw = [OH]*[H]             Vattnets protolysjämvikt
% [H] = [A] - [A]0 + [OH]                   Protonbalans (H kommer från
%                                           vatten som bildar OH och HA som bildar A, eller omvänt)

%Redigera dessa värden för olika system:
nHA0=1e-4;                      %Ämnesmängd tillsatt syra HA (i mol).
nA0=0;                          %Ämnesmängd tillsatt konjugerad bas A (i mol).
V=1;                            %Volym av vattenlösningen (i liter).
pKa=9.886;                      %Detta är pKa för fenol.


%Gissningsvärden från argumentet x1
HA=x1(1);                       %Koncentrationen av syran HA i mol/liter vid jämvikt
A=x1(2);                        %Koncentrationen av den konjugerade basen A i mol/liter vid jämvikt
OH=x1(3);                       %Hydroxidjonkoncentrationen i mol/liter vid jämvikt
H=x1(4);                        %Vätejonkoncentrationen i mol/liter vid jämvikt


pKw=14;         
Ka=10^(-pKa);
Kw=10^(-pKw);                   %Vattnets jonprodukt vid 25C

%Residualen b=f(x):
b=       [  Ka*HA-A*H;          %Syrans protolysjämvikt
            nHA0+nA0-V*(HA+A);  %Balans av HA 
            Kw-OH*H;            %Vattnets protolysjämvikt
            H+nA0/V-(A+OH)];    %Balans av H

%Jacobimatrisen B=Df(x):       
B=       [  Ka  -H  0   -A;     
            -V  -V  0   0;
            0   0   -H  -OH;
            0   -1  -1  1];
        
h=B\b;                 %Lös den linjäriserade ekvationen Bh=b
x2=x1-h;               %Uppdatera

%Om pH är bestämt och vätejonskoncentrationen ges av x1(4)
%stryks ekvation 3 och 4 genom att två föjande rader läggs till programmet:
%x2=x1;
%x2(1:2)=x1(1:2)-B(1:2,1:2)\b(1:2);
pH=-log10(x2(4));
