A method for magnetic resonance imaging, MRI, of spin systems using an MRI device, comprising - applying a consecutive series of MRI sequences (P n,P n+1,P n+2) to a target volume in the MRI device, wherein each MRI sequence (Pn) comprises an ordered combination of radio frequency excitation pulses and/or magnetic field gradient pulses configured according to experimental settings (T R, α, β) for manipulating a magnetization state (Mn) of the spin systems to be imaged in the target volume, wherein the MRI sequences are applied during a transient phase with a repetition time between subsequent sequences shorter than a relaxation time of the spin systems to be imaged; - measuring a discrete sequence of transient response signals (S n, S n+1,S n+2), wherein a respective response signal (Sn) is acquired during a respective acquisition windows of a respective MRI sequence (Sn), wherein the sequence of transient response signals represents sequential evolution of the magnetization state (Mn) resulting from a combination of the manipulating of the spin systems by the MRI sequences as well as intrinsic properties (PD, T1,T2) of the spin systems to be imaged, wherein the transient response signals are measured in respective acquisition windows per MRI sequence during the transient phase before a steady state response signal to the consecutive series of MRI sequences develops; - fitting the discrete sequence of transient response signals to a fit function that is based on an analytically modeled evolution of the magnetization state, wherein the magnetization state is modeled as a three dimensional magnetization vector, wherein an effect of an MRI sequence on the magnetization is modeled using matrices as operators on the magnetization vector, wherein consecutive application of MRI sequences is modelled as repeated multiplication by an operator matrix, wherein the repeated matrix multiplication is calculated using eigenvalues of the operator matrix, wherein the modelled fit function is based on a projection of the magnetization vector in a measurement direction of the measured responses, wherein a shape of the fit function is determined according to the analytically modelled evolution based at least on the experimental settings (T R, α, β) as well as intrinsic parameters (r,λ 3,φ,δ) to be fitted, wherein the intrinsic parameters are related to the intrinsic properties (PD,T1,T2) of the spin systems to be imaged; and - determining an image of the target volume based on one or more fitted values of the variable intrinsic parameters which match the shape of the fit function to the discrete sequence of transient response signals. - applying a consecutive series of MRI sequences (P n,P n+1,P n+2) to a target volume in the MRI device, wherein each MRI sequence (Pn) comprises an ordered combination of radio frequency excitation pulses and/or magnetic field gradient pulses configured according to experimental settings (T R, α, β) for manipulating a magnetization state (Mn) of the spin systems to be imaged in the target volume, wherein the MRI sequences are applied during a transient phase with a repetition time between subsequent sequences shorter than a relaxation time of the spin systems to be imaged; - measuring a discrete sequence of transient response signals (S n, S n+1,S n+2), wherein a respective response signal (Sn) is acquired during a respective acquisition windows of a respective MRI sequence (Sn), wherein the sequence of transient response signals represents sequential evolution of the magnetization state (Mn) resulting from a combination of the manipulating of the spin systems by the MRI sequences as well as intrinsic properties (PD, T1,T2) of the spin systems to be imaged, wherein the transient response signals are measured in respective acquisition windows per MRI sequence during the transient phase before a steady state response signal to the consecutive series of MRI sequences develops; - fitting the discrete sequence of transient response signals to a fit function that is based on an analytically modeled evolution of the magnetization state, wherein the magnetization state is modeled as a three dimensional magnetization vector, wherein an effect of an MRI sequence on the magnetization is modeled using matrices as operators on the magnetization vector, wherein consecutive application of MRI sequences is modelled as repeated multiplication by an operator matrix, wherein the repeated matrix multiplication is calculated using eigenvalues of the operator matrix, wherein the modelled fit function is based on a projection of the magnetization vector in a measurement direction of the measured responses, wherein a shape of the fit function is determined according to the analytically modelled evolution based at least on the experimental settings (T R, α, β) as well as intrinsic parameters (r,λ 3,φ,δ) to be fitted, wherein the intrinsic parameters are related to the intrinsic properties (PD,T1,T2) of the spin systems to be imaged; and - determining an image of the target volume based on one or more fitted values of the variable intrinsic parameters which match the shape of the fit function to the discrete sequence of transient response signals.