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 encontrar 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 propõe-se um novo método para a simulação de uma classe de sistemas vibratórios estocásticos. Trata-se de 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 da exponencial de uma matriz por cada passo de integração para sua implementação. Demonstra-se que o método 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.