牛顿型方法是数值求解非线性方程组的标准方法。为了确保他们能够处理这个问题，我尝试在 julia 中使用NLsolve.jl解决您的问题，它似乎工作得很好：

using NLsolve
using LineSearches
using ForwardDiff

function solve(d::Real, μ::Real, β::Real, T::Integer)
    @assert 0 < d < 1
    @assert 0 < μ + β < 1
    n = T + 1
    rhs(p₋, p, p₊) = -p + d + (1-d)*p*((1-μ-β)*p+μ*p₊+β*p₋)
    function f!(p, val)
        val[1] = p[1]
        for k = 2:T
            val[k] = rhs(p[k-1], p[k], p[k+1])
        end
        val[n] = p[n] - 1
    end
    function g!(p, J)
        J[1,1] = J[n,n] = 1
        for k = 2:T
            J[k,k-1:k+1] = ForwardDiff.gradient(x -> rhs(x...), p[k-1:k+1])
        end
    end
    Jₚ = spzeros(n, n)
    Jₚ[1,1] = Jₚ[n,n] = 1
    for k = 2:T; Jₚ[k,k-1:k+1] = 1; end
    p₀ = ones(n) * min(1, d/(1-d))
    df = DifferentiableGivenSparseMultivariateFunction(f!, g!, Jₚ)
    @time sol = nlsolve(df, p₀, method=:newton, linesearch! = LineSearches.strongwolfe!, show_trace = true)
    g!(sol.zero, Jₚ)
    @printf "Jacobian condition number: %.2e\n" cond(Array(Jₚ), Inf)

    sol
end

我注意到这个案子$\mu+\beta\approx 1$这里似乎是最困难的情况，解决方案取决于一个好的初始猜测，因为方程是刚性的，并且在$0$和$T$. 所以这不是一个完整的答案，它只适用于非极端情况。

一般来说，不能指望这样一个问题有一个封闭形式的解决方案。此外，Groebner 基础在相对较小的系统上是可行的，而不是像这个系统。确实，这里的雅可比行列式是稀疏的，但这还不够。简单地拥有一个稀疏的非线性方程组并不一定意味着会有一种简单有效的方法来解决它们。

此外，正如关于 MO 的这个问题所提到的，二次方程的多元系统通常是一种非常困难的问题，至少是 NP-hard 并且可能更难。

除了基里尔的出色回答外，我认为工作提到的一件事是斯梅尔的$\alpha$-理论。应用牛顿方法的难点之一是，如果您没有做出非常好的初始猜测，迭代可能会发散。为了解决这个问题，通常使用阻尼牛顿迭代，但要选择阻尼参数，您所拥有的只是启发式方法。

Smale 的理论告诉您，当函数是解析函数时，牛顿法的初始猜测何时位于非线性方程的某个根的二次收敛区域内的充分条件。充分条件要求您对所讨论函数的泰勒系数有统一的界限。对于多项式，条件特别容易评估，因为泰勒系数在多项式的次数之后为 0，在您的情况 2 中。此外，您正在解决的系统的稀疏性使得评估应用 Smale 理论的条件更加简单。

我的老同事有一些幻灯片，我认为这些幻灯片非常适合作为该主题的参考，他也有一些Python示例代码。