The Thomas Method Algorithm for solving three diagonal system in Multi-step Wide-Angle (Pade(2,2)) Beam Propagation Method

Question

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}$$

1. 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