I'm trying to simulate in Matlab a 2D gaussian beam propagation in free propagation region using multi-step Wide-Angle BPM (Pade(2, 2)) with Crank-Nicholson scheme. According to book Pade polynomials coefficients are given by:

$$\xi_1 = \frac{3}{4}-\frac{iK\Delta z(1-\alpha)}{2}$$
$$\xi_2 = \frac{1}{16}-\frac{iK\Delta z(1-\alpha)}{4}$$
$$\chi_1 = \frac{3}{4}+\frac{iK\Delta z\alpha}{2}$$
$$\chi_2 = \frac{1}{16}+\frac{iK\Delta z\alpha}{4}$$

where $\alpha$ is backward/forward component ratio and $\alpha=0.5$.

Roots of the numerator and denominator polynomials can be obtained:

$$\alpha_1=\frac{1}{2}\left[\xi_1-\sqrt{\xi_1^2-4\xi_2}\right]$$
$$\alpha_2=\frac{1}{2}\left[\xi_1+\sqrt{\xi_1^2-4\xi_2}\right]$$
$$\beta_1=\frac{1}{2}\left[\chi_1-\sqrt{\chi_1^2-4\chi_2}\right]$$
$$\beta_2=\frac{1}{2}\left[\chi_1+\sqrt{\chi_1^2-4\chi_2}\right]$$

The propagation is split into two steps:

$$u^{m+\frac{1}{2}}=\frac{(1+\alpha_1\mathcal{P})}{(1+\beta_1\mathcal{P})}u^m$$
$$u^{m+1}=\frac{(1+\alpha_2\mathcal{P})}{(1+\beta_2\mathcal{P})}u^{m+\frac{1}{2}}$$

which yields finite-difference equation

$$u_j^{m+\frac{1}{2}}+\frac{\beta_1}{K^2}\frac{u_{j-1}^{m+\frac{1}{2}}-2u_j^{m+\frac{1}{2}}+u_{j+1}^{m+\frac{1}{2}}}{\Delta x^2}+\frac{\beta_1}{K^2}(k^2-K^2)u_j^{m+\frac{1}{2}}=u_j^m+\frac{\alpha_1}{K^2}\frac{u_{j-1}^m-2u_j^m+u_{j+1}^{m}}{\Delta x^2}+\frac{\alpha_1}{K^2}(k^2-K^2)u_j^m$$

where $k = k_0n$ is position-dependent wavenumber and $K=k_0n_0$ is the reference wavenumber. This equation forms a tridiagonal system, from where its coefficient $a_j, b_j, c_j, r_j$ are given by:

$$a_j=\frac{\beta_1}{K^2\Delta x^2}$$
$$b_j=1-\frac{2\beta_1}{K^2\Delta x^2}+\frac{\beta_1}{K^2}(k^2-K^2)$$
$$c_j=\frac{\beta_1}{K^2\Delta x^2}$$
and $r_j$ is the right side of equation.

This is tridiagonal system I need to solve. To do this I use Thomas Method Algorithm. This method says:

1. Set $\beta=b_1, u_1=\frac{r_1}{\beta}$

2. Evaluate for j=2 to n:

$$\gamma_j=\frac{c_{j-1}}{\beta}$$
$$\beta=b_j-a_j\gamma_j$$
$$u_j=\frac{r_j-a_ju_{j-1}}{\beta}$$

- Find for j=1 to n-1

$$k=n-j$$
$$u_k=u_k-\gamma_{k+1}u_{k+1}$$

Now I have problem, because I can't really figure out how to correctly implement Thomas algorithm. Since in this situation $a_j, b_j, c_j = const$ only $r_j$ changes but I don't understand 1st point in this method. This is my Matlab code:

```
function u_m_total = WideAnglePade2(k0, step_z, step_x, lateral_x, ...
input_field, z_steps, n)
%% COEFFICIENTS
alpha = 0.5;
xi_1 = 0.75 - (1j * k0 * step_z * (1 - alpha))/2;
xi_2 = 1/16 - (1j * k0 * step_z * (1 - alpha))/4;
chi_1 = 0.75 + (1j * k0 * step_z * alpha)/2;
chi_2 = 1/16 + (1j * k0 * step_z * alpha)/4;
alfa1 = 0.5 * (xi_1 - sqrt(xi_1^2 - 4 * xi_2));
alfa2 = 0.5 * (xi_1 + sqrt(xi_1^2 - 4 * xi_2));
beta1 = 0.5 * (chi_1 - sqrt(chi_1^2 - 4 * chi_2));
beta2 = 0.5 * (chi_1 + sqrt(chi_1^2 - 4 * chi_2));
% caching coefficients to speed up computation
coeff_a1 = beta1 / (k0^2 * step_x^2);
coeff_b1 = 1 - 2 * beta1 / (k0^2 * step_x^2) + beta1 / k0^2 * ((k0 * n)^2 - k0^2);
coeff_c1 = beta1 / (k0^2 * step_x^2);
coeff_a2 = beta2 / (k0^2 * step_x^2);
coeff_b2 = 1 - 2 * beta2 / (k0^2 * step_x^2) + beta2 / k0^2 * ((k0 * n)^2 - k0^2);
coeff_c2 = beta2 / (k0^2 * step_x^2);
beta_1 = coeff_b1;
beta_2 = coeff_b2;
% creating vectors holding values for m_half_z and m_z
% I have problem in here, I don't know what to assign and where
u_m_half = zeros(z_steps, lateral_x);
u_m_total(1, :) = input_field;
u_m_total(1, 1) = input_field(1) / beta2;
for m=1:z_steps-1
disp(m);
%% STEP 1 PROPAGATION ΔZ/2
gamma = coeff_c1 / beta_1;
beta = coeff_b1 - coeff_a1 * gamma;
for j=2:lateral_x-1
r = u_m_total(m, j) + alfa1/k0^2 * (u_m_total(m, j-1) - ...
2*u_m_total(m,j) + u_m_total(m, j+1))/step_x^2 + alfa1/k0^2 *...
((k0 * n)^2 - k0^2) * u_m_total(m, j);
u_m_half(m, j) = (r - coeff_a1 * u_m_total(m, j-1))/beta;
end
for j=1:lateral_x-1
k=lateral_x-j;
u_m_half(m, k) = u_m_half(m, k) - gamma * u_m_half(m, k+1);
end
%% STEP 2 PROPAGATION ΔZ
gamma = coeff_c2 / beta_2;
beta = coeff_b2 - coeff_a2 * gamma;
for j=2:lateral_x-1
r = u_m_half(m, j) + alfa2/k0^2 * (u_m_half(m, j-1) - ...
2*u_m_half(m,j) + u_m_half(m, j+1))/step_x^2 + alfa2/k0^2 ...
* ((k0 * n)^2 - k0^2) * u_m_half(m, j);
u_m_total(m+1, j) = (r - coeff_a2 * u_m_half(m, j-1))/beta;
end
for j=1:lateral_x-1
k=lateral_x-j;
u_m_total(m+1, k) = u_m_total(m+1, k) - gamma * u_m_total(m+1, k+1);
end
end
end
```

But instead of simulating dispersion and broadening of beam, its amplitude goes towards +Inf.

This post imported from StackExchange Physics at 2017-01-11 10:31 (UTC), posted by SE-user Colonder