# Solving the PDE of a signal in an optical fiber

Disclaimer Please let me know if this is not the right place to ask this question.
Disclaimer Part 2 I have no background in signals. My focus is in mathematics.

In my programming class, we are studying the pulse propagation in an optical fiber. The fiber can be modeled as a line marked by a one-dimensional coordinate $$z$$. The pulse at the fiber input ($$z=0$$) is given as a time-dependent field envelope $$A(0,t)$$. I need to derive the magnitude of the field envelope at any location and at any time $$A(z,t)$$. We are told to use numerical methods (i.e. Fourier transform) to solve the partial differential equation

$$\dfrac{\partial A}{\partial z} = \dfrac{\alpha}{2} A - \beta_1\dfrac{\partial A}{\partial t} -j\dfrac{\beta_2}{2}\dfrac{\partial^2 A}{\partial t^2}$$

Here, the first term on the right-hand side models the signal attenuation due to scattering, the second term models the group velocity, and the third term models chromatic dispersion (i.e. different colors travel at different speeds in the fiber, resulting in pulse broadening). My teacher has said to use a Gaussian pulse and to discretize the pulse over the range $$\left[-\frac{t}{2}, \frac{t}{2}\right]$$ such that $$t >> \sigma$$ $$\left(t=10\sigma\right)$$.

I have been able to show mathematically that the solution should be of the form $$A(z,t) = \text{ifft}\left[\tilde{A}(0,\omega)\exp\left(\left(\frac{\alpha}{2} - j\omega\beta_1 + \frac{j\omega^2}{2}\beta_2\right)z\right)\right]$$ (where $$\tilde{A}(0,\omega)$$ is the Fourier transform of $$A(0,t)$$) but am having trouble actually implementing these results in MATLAB.

What I have so far:

sigma = 0.1; % standard deviation of the pulse
N = 256;  % number of points sampled
t = linspace(-10*sigma, 10*sigma, N); % time vector
A0t = (0.25/sqrt(2*pi*sigma^2))*exp(-(t.^2)/(2*sigma^2)); % descritized gaussian pulse

Fs = 1000; %sampling frequency


(I don't know what this number should be)

A0w = fft(A0t);
A0w = abs(A0w(1:round(N/2)));

w = (0:(N-1))*(Fs/N);  %frequency vector
w = w(w < Fs/2);


(I don't know why the last three lines are necessary. I saw it here.)

alpha = input(...);
beta1 = input(...);
beta2 = input(...);

lambda = ((alpha/2) - (1j*w*beta1) + (1j*(w.^2)*beta2/2)); %eigenvalue of ODE in frequency domain

end = input(...);

z = linspace(0,end,length(A0w))'; % spatial vector

Azw = A0w.*exp(lambda.*z); % solution in the frequency domain

Azt = abs(ifft2(Azw)); % solution in the time domain


Using alpha = 10, beta1 = 50, beta2 = 18, and end = 10 results in the surface plot below, which I do not believe is correct. ($$t$$ and $$z$$ are plotted against their indices rather than the actual values because $$t$$ is not the same size as $$A(z,t)$$.)

If somebody could please explain to me why I need to do the things that I am doing in my code (the italic comments) as well as what I should be doing otherwise to get a correct output, or if what I am getting is correct, it would be much appreciated. 