n=4;
x1=zeros(1,n);
x2=(n+1)*ones(1,n);
x3=[1:n];
x=[x1,x2,x3];
y1=[1:n]/n;
y2=y1;
y0=[1:n]/(n+1)
y3=1+1*(y0).*(1-y0);
y=[y1,y2,y3];
figure;
%scatter(x,y)
tri3=zeros(2*n,3);
for k=1:n-1
tri3(1+2*(k-1),:)=[k,k+1,n+k];
tri3(2+2*(k-1),:)=[k+1,n+k,n+k+1];
end % for k
for k=1:n-1
tri3(k+2*n-2,:)=[n,2*n+k,2*n+k+1]
end % for k
tri3(3*n-2,:)=[n,2*n,3*n];
plot_tri(x,y,tri3)
M=make_adj_matrix(tri3);
B=find_boundary(tri3,M);
plot_tri_boundary(x,y,B);
radii=find_radii2(x,y,tri3);
plot_central_regions2(x,y,tri3,radii);
propagate2(x,y,tri3,radii,B)
%draw_prop(x,y,tri3)