t0=linspace(0,20,21);
t=1790:10:1900;
x=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76.0,92.0,106.5,123.2,131.7,150.7,179.3,204.0,226.5,251.4,281.4];
a0=[392.0886 0.2557];
f=@(a,t) a(1)./(1+a(1)/x(1)-1)*exp(-a(2)*t0);
a=lsqcurvefit(f,a0,t,x)
y=f(a,t);
plot(t0,x,'+',t0,y,'b-')
