function [b0,b1,rsq,StdErrEst] = weighted_linfit(x,y,dy)
n=length(x);
m=1;
V=ones(n,2); %initalize V matrix
V(:,2)=x;
weights=1./dy;
W=diag(weights);
Wy=W*y;
WV=W*V;
A=WV'*WV;
b=WV'*Wy;
w=A\b; %solve Normal Equations
b0=w(1);
b1=w(2);
fit=V*w;
plot(x,fit)
hold on
errorbar(x,y,dy,'*');
hold off
%standard error of estimataion and coefficient of determination
r=fit-y;
sse=r'*r;
StdErrEst=sqrt(sse/(n-m-1));
ybar=sum(y)/n;
var=y-ybar;
sst=var'*var;
ssr=sst-sse;
rsq=ssr/sst;
end