早速コード


function newmarkβ(input, Δt, ω, h, β=0.25, ε=0.000001)
  output = zeros(input)

  u̇ₓ = u̇₁ = zeros(input[:, 1])
  uₓ = u₁ = output[:, 1] = zeros(input[:, 1])
  üₓ = ü₁ = -input[:, 1] - 2 * h * ω * u̇₁ - ω ^ 2 * u₁

  for x=1:length(input)-1
    ǖₓ₊₁ = üₓ
    while true
      u̇ₓ₊₁ = u̇ₓ + 0.5 * (üₓ + ǖₓ₊₁) * Δt
      uₓ₊₁ = uₓ + u̇ₓ * Δt + (0.5 - β) * üₓ * Δt ^ 2 + β * ǖₓ₊₁ * Δt ^ 2

      üₓ₊₁ = -input[:, x+1] - 2 * h * ω * u̇ₓ₊₁ - ω ^ 2 * uₓ₊₁

      if norm(ǖₓ₊₁ - üₓ₊₁) < ε
        üₓ = üₓ₊₁
        u̇ₓ = u̇ₓ₊₁
        uₓ = uₓ₊₁
        output[:, x+1] = uₓ₊₁
        break
      end
      ǖₓ₊₁ = üₓ₊₁
    end
  end
  return output
end

t = transpose(linspace(0.0, 100, 1000))
Δt = t[2]
input = sin(t)
res=newmarkβ(input, Δt, 1.0, 0.1)
using Gadfly
plot(x=t, y=res)

擬似コードに見せかけたJuliaで実行可能なコード
ユニコードプログラミングをネタ的にやってみた
普通に変位応答を戻り値で返しているが、実際にはobserverを受け取る実装にしたほうが良いだろう。
ここでは線形の範囲（多自由度にも一応対応しているはず）で解いているので全くNewmark-βの意味が無い。
非線形応答をやるなら剛性と減衰を関数か何かで渡す必要がある。コルーチンがあれば楽かな多分

自分用：Newmark-β法

早速コード