function [x,h] = belanger(Q,C,L,B,ib,xb,hb,N);

% Set constant
g         = 9.81;

% Compute he and hg
hg3       = (Q.^2./B.^2./g)
he3       = (Q.^2./B.^2./C.^2./ib)
hg        = hg3.^(1/3)

% Compute cf
cf        = g./(C.^2);

% Define grid
dx        = L/N
x         = [xb-L:dx:xb];
sx        = length(x);

% Compute surface level
h(sx)     = hb;
for i = (sx-1):-1:1;
    
    if ib == 0;
        %rhs = -(cf*hg.^3)./(h(i+1).^3-hg3);
        rhs = -(cf*hg.^3)./(h(i+1).^3-hg3);
    else
        rhs = ib.*(h(i+1).^3-he3)./(h(i+1).^3-hg3);
    end
    h(i)  = h(i+1) - dx.*rhs;
    
end
