function b=bez_spline(u,b0,v0,vn);
  k=length(u)-1; %stevilo segmentov
  d=zeros(1,k);
  for i=1:k
    d(i)=u(i+1)-u(i); %delta_i
  end
  m=2*k-2;
  A=zeros(m);
  b=zeros(2,3*k+1);
  b(:,[1:3:3*k+1])=b0;
  b(:,2)=b(:,1)+d(1)*v0;
  b(:,3*k)=b0(:,end)-d(k)*vn
  A(1:2,:)=[[d(2),d(1),zeros(1,m-2)];[2*d(2)^2,-2*d(1)^2,d(1)^2,...
  zeros(1,m-3)]];
  A(m-1:m,:)=[[zeros(1,m-2),d(k),d(k-1)];[zeros(1,m-3),...
  -d(k)^2,2*d(k)^2,-2*d(k-1)^2]];
  for i=3:2:m-2
    A([i,i+1],[i-1:i+2])=[[0,d((i+3)/2),d((i+3)/2-1),0];...
    [-d((i+3)/2)^2,2*d((i+3)/2)^2,...
    -2*d((i+3)/2-1)^2,d((i+3)/2-1)^2]];
  end;
  A
  y=zeros(m,2);
  y(2,:)=(d(2)^2-d(1)^2)*b(:,4)'+d(2)^2*b(:,2)';
  y(m,:)=(d(k)^2-d(k-1)^2)*b(:,end-3)'-d(k-1)^2*b(:,end-1)';
  
  for i=4:2:m-2
    y(i,:)=(d((i+2)/2)^2-d((i+2)/2-1)^2)*b(:,3*(i)/2+1)';
  end;
  b
  for i=1:2:m-1
    y(i,:)=(d((i+1)/2+1)+d((i+1)/2))*b(:,3*(i+1)/2+1)';
  end
  y
  bb=A\y;
  b(:,3)=bb(1,:)';
  b(:,end-2)=bb(end,:)';
  for i=2:k-1
    b(:,3*i-1:3*i)=bb([2*i-2:2*i-1],:)';
  end;
return
