def solve(self, u, b, solvers1D, naxes):
        if u is not b:
            u[:] = b

        s = [0]*u.ndim
        s[naxes] = slice(None)
        paxes = np.setxor1d(range(u.ndim), naxes)
        if u.ndim == 2:
            for i, sol in enumerate(solvers1D):
                s[paxes[0]] = i
                s0 = tuple(s)
                sol.inner_solve(u[s0], sol._inner_arg)

        elif u.ndim == 3:
            for i, m in enumerate(solvers1D):
                s[paxes[0]] = i
                for j, sol in enumerate(m):
                    s[paxes[1]] = j
                    s0 = tuple(s)
                    sol.inner_solve(u[s0], sol._inner_arg)