function O = oxygen_solve(H,vb,delta_z,O_bound,D_o,delta_ox,ex)

    A = zeros(H/delta_z+1,H/delta_z+1);
    B = [zeros(H/delta_z,1);O_bound];
    for i = 1:(H/delta_z+1)
        if i == 1
            A(i,i) = delta_z^2*delta_ox/D_o*ex(i) + 2;
            A(i,i+1) = -2;
        elseif i == H/delta_z+1
            A(i,i) = 1;
        else
            A(i,i) = delta_z^2*delta_ox/D_o*ex(i) + 2;
            A(i,i-1) = -1;
            A(i,i+1) = -1;
        end
    end
    O = linsolve(A,B);
    O = [O;repmat(O_bound,vb/delta_z,1)];

end