x = newton('vasafunk',[0 1e-7 0 0 0 0 0 0 1e-7 ]',1e-15)

pH = 5:13;
c1 = [0.1 0.1 0.1 0.09999999999993 0.09999999999272 0.09999999893183 0.09999955228507 0.09961580108601 0.02085605080772];
c2 = [0.09999999893009 0.09999998313902 0.09999921071079 0.099929746221630.09321266629094 0.00856041152907 0.00002233055330 0.00000002592818 0.00000000002635];

ph = linspace(5,13,100);
C1 = spline(pH,c1,ph);
C2 = spline(pH,c2,ph);

plot(ph,C1,pH,c1,'+',ph,C2,pH,c2,'+')
axis([5 13 0 0.15 ])
xlabel('pH')
legend('Fe(EDDHMA)^-','Fe(EDTA)^-')
ylabel('M')
grid