For testing I would like to solve the equation$$\partial_tA=\nabla^2_xA,\,A_0=\sin(x)$$In order to solve it, I use an implicit solver$$\left(1-\frac{\nabla^2_x}{2}dt\right)A_{n+1} = \left(1+\frac{\nabla^2_x}{2}dt\right)A_n$$or the Runge-Kutta-method in two different orders. For solving the whole equation I use the following matlab script:clear all;warning('off', 'all');x_max = 1;x_num = 1000;a0_val = 1;b0_val = 1;dx = x_max/20000;steps = 1/dx;x_vec = linspace(0, x_max, x_num).';A0_vec = sin(x_vec);%linspace(1, x_num, x_num).';ones_vec = ones(x_num, 1);dh = (abs(x_vec(2)-x_vec(1)));nabla_matrix = spdiags([-1*ones_vec 16*ones_vec -30*ones_vec 16*ones_vec -1*ones_vec], -2:2, x_num, x_num);nabla_matrix(1, 1) = 35;nabla_matrix(1, 2) = -104;nabla_matrix(1, 3) = 114;nabla_matrix(1, 4) = -56;nabla_matrix(1, 5) = 11;nabla_matrix(2, 1) = 11;nabla_matrix(2, 2) = -20;nabla_matrix(2, 3) = 6;nabla_matrix(2, 4) = 4;nabla_matrix(2, 5) = -1;nabla_matrix(x_num, x_num) = 35;nabla_matrix(x_num, x_num-1) = -104;nabla_matrix(x_num, x_num-2) = 114;nabla_matrix(x_num, x_num-3) = -56;nabla_matrix(x_num, x_num-4) = 11;nabla_matrix(x_num-1, x_num) = 11;nabla_matrix(x_num-1, x_num-1) = -20;nabla_matrix(x_num-1, x_num-2) = 6;nabla_matrix(x_num-1, x_num-3) = 4;nabla_matrix(x_num-1, x_num-4) = -1;nabla_matrix = nabla_matrix / (12*(dh*dh));A_equation = @(A_vec) ((nabla_matrix*A_vec.'));A_matrix_single = speye(x_num)-dx/2.*nabla_matrix;B_matrix_single = speye(x_num)+dx/2.*nabla_matrix;A_matrix = zeros(steps, x_num);A_matrix_NL_rk4 = zeros(steps, x_num);A_matrix_euler = zeros(steps, x_num);A_matrix_rk4 = zeros(steps, x_num);A_matrix_rk42 = zeros(steps, x_num);A_matrix(1,:) = A0_vec;A_matrix_NL_rk4(1,:) = A0_vec;A_matrix_euler(1,:) = A0_vec;A_matrix_rk4(1,:) = A0_vec;A_matrix_rk42(1,:) = A0_vec;%Do stepsfor i = 2:steps    A_k1 = A_equation(A_matrix_rk4(i-1,:));    A_k2 = A_equation(A_matrix_rk4(i-1,:)+A_k1.'*dx/2);    A_k3 = A_equation(A_matrix_rk4(i-1,:)-A_k1.'*dx+2*dx*A_k2.');    A_matrix_rk4(i,:) = A_matrix_rk4(i-1,:)+dx*(A_k1/6+2/3*A_k2+A_k3/6).';endfor i = 2:steps    A_k1 = A_equation(A_matrix_rk42(i-1,:));    A_k2 = A_equation(A_matrix_rk42(i-1,:)+A_k1.'*dx/2);    A_k3 = A_equation(A_matrix_rk42(i-1,:)+0.5*dx*A_k2.');    A_k4 = A_equation(A_matrix_rk42(i-1,:)+dx*A_k3.');    A_matrix_rk42(i,:) = A_matrix_rk42(i-1,:)+dx*(A_k1/6+1/3*A_k2+A_k3/3+A_k4/6).';endfor i = 2:steps    A_k1 = A_equation(A_matrix_NL_rk4(i-1,:));    A_k2 = A_equation(A_matrix_NL_rk4(i-1,:)+A_k1.'*dx/2);    A_k3 = A_equation(A_matrix_NL_rk4(i-1,:)+0.5*dx*A_k2.');    A_k4 = A_equation(A_matrix_NL_rk4(i-1,:)+dx*A_k3.');    A_matrix_NL_rk4(i,:) = speye(x_num)\(speye(x_num)*A_matrix_NL_rk4(i-1,:).'+speye(x_num)*(dx*(A_k1/6+1/3*A_k2+A_k3/3+A_k4/6)));endfor i = 2:steps    A_matrix(i,:) = A_matrix_single\(B_matrix_single*A_matrix(i-1,:).');endfigure;imagesc(A_matrix);title('Forward step')figure;imagesc(A_matrix_rk4);title('RK4 step')figure;imagesc(A_matrix_rk42);title('RK42 step')figure;imagesc(A_matrix_NL_rk4);title('RK4_NL step')but the results are way of my expectations. I get a more or less correct result when using the implicit approach, but the latter three all explode.Why is that? Is my general approach flawed?
Comparison of different methods for solving $\partial_tA=\nabla^2_xA$
matlab;pde