(*
x=y=Table[i,{i,1,10}];
sx=sy=Table[i/10,{i,1,10}];
*)

x={1.5,3,4,5,6,7,8,9,10};
sx={1,0.5,0.5,0.5,0.5,0.5,0.5,1.5,1.5};
y={0.5,1,2,3,4,5.5,7,8.5,10};
sy={1,1,1,0.5,0.5,0.5,1,1,1.5};

r=0.3; (* correlation between errors in x and y, if it exists *)
n=Length[x];
wx=wy=W=u=v=β=ξ=η=Table[0,{n}];

a1=1.0; (* initial approx for slope *)

For[i=1,i≤n,i++,{
wx[[i]]=1/((sx[[i]])^2);
wy[[i]]=1/((sy[[i]])^2);
}];

(* iteration to get the slope a1 *)
For[k=1,k≤10,k++,{
For[i=1,i≤n,i++,{
W[[i]]=wx[[i]]wy[[i]]/(wx[[i]]+(a1^2)wy[[i]]-2a1 r Sqrt[wx[[i]]wy[[i]]]);
}];
xave=Sum[W[[i]]x[[i]],{i,1,n}]/Sum[W[[i]],{i,1,n}];
yave=Sum[W[[i]]y[[i]],{i,1,n}]/Sum[W[[i]],{i,1,n}];
For[i=1,i≤n,i++,{
u[[i]]=x[[i]]-xave;
v[[i]]=y[[i]]-yave;
β[[i]]=W[[i]](u[[i]]/wy[[i]]+a1 v[[i]]/wx[[i]]-(a1 u[[i]]+v[[i]])r/Sqrt[wx[[i]]wy[[i]]]);
}];
a1=Sum[W[[i]]β[[i]]v[[i]],{i,1,n}]/Sum[W[[i]]β[[i]]u[[i]],{i,1,n}];
Print[k," ",a1]
}];

(* solution for intercept *)
a0=yave-a1 xave;

For[i=1,i≤n,i++,ξ[[i]]=xave+β[[i]];];
xiave=Sum[W[[i]]ξ[[i]],{i,1,n}]/Sum[W[[i]],{i,1,n}];
For[i=1,i≤n,i++,η[[i]]=ξ[[i]]-xiave;];
sig2a1=1/Sum[W[[i]]((η[[i]])^2),{i,1,n}];
sig2a0=1/Sum[W[[i]],{i,1,n}]+(xiave^2)sig2a1;

Print[a0," ",a1," ",Sqrt[sig2a0]," ",Sqrt[sig2a1]]