% Forward Euler
lambda=complex(-1,0);
% stability limit for (-1,0): dt<2
% stability limit for (-1,5): dt<1/13~0.077
dt=2.1;
G=1+lambda*dt;
xe(1)=1;
xi(1)=1;
t(1)=0;
tend=20;
Nsteps=round(tend/dt);
%Forward Euler
for i=2:Nsteps
xe(i)=G*xe(i-1);
% xe(i)=exp(lambda*dt)*xe(i-1);
t(i)=t(i-1)+dt;
end
% plot(t,real(xe),'b',t,real(x),'m')
hold on;
% pause
% Backward Euler
G=1/(1-lambda*dt);
for i=2:Nsteps;
xi(i)=G*xi(i-1);
% xe(i)=exp(lambda*dt)*xe(i-1);
% t(i)=t(i-1)+dt;
end
% plot(t,real(x),'r')
%
% generate exact solution with small time steps
tlist=0:tend/1000:tend;
xex=real(exp(lambda*tlist));
plot(t,real(xe),'-o',t,real(xi),'-x',tlist,xex,'-');
xlabel('t');
ylabel('x(t)');
legend('EE','IE','Exact');
data1 = [t', xe', xi'];
data2 = [tlist', xex'];
save "TimeIntegrationExample.dat" data1 data2