A method for processing seismic data comprising: a) using seismic sources to generate seismic waves that propagate into the earth; b) obtaining seismic data by detecting with a seismic receiver the signal from reflected seismic waves; and c) processing, in a computer, the seismic data to generate a prediction of internal multiples at the receiver location and outputting predicted internal multiples for attenuation; characterized by: d) generating internal multiples by solving a series of partial differential wave equations, wherein a first partial differential wave equation describes propagation of a seismic wave going from a first reflector to a second reflector, wherein a second partial differential wave equation describes propagation of a seismic wave going from the second reflector to a third reflector, and wherein a third partial differential wave equation describes propagation of a seismic wave going from the third reflector to a seismic receiver; wherein the first partial differential equation is ∂ ∂ z + iq p ^ 1 k ω z k s = ∫ b ^ 1 k; − k → ′; z p ^ 0 k → ′; ω; z; k s d k → ′, wherein z is depth between z 1 and z 2, p̂ 0 is Fourier transform of p 0 which is Fourier conjugate to field point going from source to the reflections above depth z 1, p̂ 1 is Fourier transform of p 1 which is Fourier conjugate to field point going from z 1 to z 2, b 1 is Fourier transform of b 1, b 1 is first term of inverse scattering series, k ' is wave vector directed towards observation point, and k s is Fourier conjugate of x s, ω is radial frequency, k s is Fourier conjugate of x s which is source coordinates, and q = ω c 1 − c 2 k x 2 ω 2 wherein c is velocity and k x is wavenumber. a) using seismic sources to generate seismic waves that propagate into the earth; b) obtaining seismic data by detecting with a seismic receiver the signal from reflected seismic waves; and c) processing, in a computer, the seismic data to generate a prediction of internal multiples at the receiver location and outputting predicted internal multiples for attenuation; characterized by: d) generating internal multiples by solving a series of partial differential wave equations, wherein a first partial differential wave equation describes propagation of a seismic wave going from a first reflector to a second reflector, wherein a second partial differential wave equation describes propagation of a seismic wave going from the second reflector to a third reflector, and wherein a third partial differential wave equation describes propagation of a seismic wave going from the third reflector to a seismic receiver; wherein the first partial differential equation is ∂ ∂ z + iq p ^ 1 k ω z k s = ∫ b ^ 1 k; − k → ′; z p ^ 0 k → ′; ω; z; k s d k → ′, wherein z is depth between z 1 and z 2, p̂ 0 is Fourier transform of p 0 which is Fourier conjugate to field point going from source to the reflections above depth z 1, p̂ 1 is Fourier transform of p 1 which is Fourier conjugate to field point going from z 1 to z 2, b 1 is Fourier transform of b 1, b 1 is first term of inverse scattering series, k ' is wave vector directed towards observation point, and k s is Fourier conjugate of x s, ω is radial frequency, k s is Fourier conjugate of x s which is source coordinates, and q = ω c 1 − c 2 k x 2 ω 2 wherein c is velocity and k x is wavenumber.