function [Rz,Rz2,NORMAL] = solve_natural_frame_from_lh(Rz0,Rz20,NORMAL0,mask,cf,h,mu,uh,vh)
% solve_natural_frame_from_lh
% solve natural frame from conformal factor cf and mean curvature H, plus
% Rz,Rz2,NORMAL on boundary, boundary is indicated by b

% first solve mu
% mu = solve_mu_from_lh(Rz0,Rz20,NORMAL,B2,cf,H,N,M,uh,vh);
% mu = zeros(size(cf));
% discretize the ODEs
[M,N] = size(cf);
A = discretizeODEs(cf,h,mu,mask,N,M,uh,vh);

% b is indicator of boundary
bf = false(M,N);

bx = [mask,bf,bf];
by = [bf,mask,bf];
bz = [bf,bf,mask];
b3 = [mask,mask,mask];
in3 = ~b3;

r = zeros(size(A,1),3);
x = zeros(size(A,1),1);
y = zeros(size(A,1),1);
z = zeros(size(A,1),1);
x(b3) = [Rz0(bx);Rz20(bx);NORMAL0(bx)];
y(b3) = [Rz0(by);Rz20(by);NORMAL0(by)];
z(b3) = [Rz0(bz);Rz20(bz);NORMAL0(bz)];

A2 = A(in3,in3);
r2 = r(in3,:) - A(in3,b3)*[x(b3),y(b3),z(b3)];
xyz2 = A2\r2;

x(in3) = xyz2(:,1);
y(in3) = xyz2(:,2);
z(in3) = xyz2(:,3);

x3 = reshape(x,[M*N,3]);
y3 = reshape(y,[M*N,3]);
z3 = reshape(z,[M*N,3]);

Rz = reshape([x3(:,1);y3(:,1);z3(:,1)],[M,N,3]);
Rz2 = reshape([x3(:,2);y3(:,2);z3(:,2)],[M,N,3]);
NORMAL = reshape([x3(:,3);y3(:,3);z3(:,3)],[M,N,3]);