close all; clear; clc;

%%
N = 100+1; % Number of discretization points
x = linspace(-1,1,N); h = x(2)-x(1);
[X,Y] = meshgrid(x); % x and y points on a grid

% Initial electric field (none)
%E = peaks(N);
%E = exp(-((X-0.5).^2+(Y-0.3).^2)*50);
E = zeros(N,N);

% Construct the Laplacian operator
D = spdiags(bsxfun(@times,[1 -2 1],ones(N^2,1)),[-1 0 1],N,N)./(4*h^2);
L = kron(D,speye(N)) + kron(speye(N),D); % Sum of first derivatives
lap = @(e) L*e; % Numerical Laplacian of vectorised E

% Driving frequency
f = 4;
isd = exp(-((X-0.5).^2+(Y-0.3).^2)*50);
d = @(t) isd(:)*pi*f*sin(2*pi*t*f)*(t<1/f);

v = 1; % Phase velocity
e0 = E(:); e0 = [e0; lap(e0)]; % Starting point
dedt = @(t,e) [e(N^2+1:end) + d(t); 2*v.^2.*lap(e(1:N^2))];
[t,U] = ode113(dedt,linspace(0,2.5,1e2),e0);

% Reshape back to 2d coordinates
U_square = zeros(N,N,numel(t));
for ti = 1:numel(t)
    U_square(:,:,ti) = reshape(U(ti,1:N^2),N,N);
end

%% Animate
for ti = 1:numel(t)
    imagesc(U_square(:,:,ti));
    caxis([-1 1]*0.2);
    colormap(plasma);
    axis square; grid on;
    drawnow; pause(0.01);
end
