Métodos de simulação numérica estão entre as mais usadas técnicas para calcular a resposta estocástica de sistemas de engenharia sujeitos a vibrações aleatórias. Em várias situações práticas os sistemas precisam ser integrados sobre longos intervalos de tempo e um grande número de amostras de trajetórias precisam ser geradas para atingir as quantidades de interesse. Nessas circunstâncias, os métodos convencionais usualmente apresentam comportamento explosivo e sua implementação demanda muito esforço computacional. Neste trabalho nós propomos um novo método para a simulação de uma classe de sistemas vibratórios estocásticos. Nós criamos um método simplificado de convergência fraca que precisa basicamente de uma variável aleatória discreta de três pontos e do cálculo do exponencial de uma matriz por cada passo de integração para sua implementação. Nós também provamos que ele reproduz as mesmas propriedades estatísticas que caracterizam a solução exata da equação linearizada associada ao sistema em questão. Isso faz do método uma ferramenta eficiente para simulações de prazo longo. Experimentos numéricos são implementados para ilustrar a performance prática do método.