%A script that use Euler's methods to solve a stiff problem

%define time span
t0=0; tend=3;

%set initial condition
%y0=[1; 1];
y0=[1; 1];

%timestep
dt=0.021; 

%the number of timesteps (no adaptive timesteps yet)
nt=round((tend-t0)/dt+realmin);

%initialize
clear y t y_b_euler;
y(:,1)=y0;
t(1)=t0;

lamda=[[-100 1];[0 -1/10]];

%Euler's method
II=eye(2);
for i=1:nt;
    t(i+1)=t(i)+dt;
    y(:,i+1)=(II+dt*lamda)*y(:,i);    
end

%exact solution
nt_exact=10000;
t_exact=(1:nt_exact+1)*(tend-t0)/nt_exact+t0;
clear y_exact;
for i=1:nt_exact+1;
    y_exact(:,i)=expm(lamda*t_exact(i))*y0;
end

%plotting the results
figure(1)
set(0,'DefaultAxesfontsize',20);
set(0,'Defaultlinelinewidth',1.2);

hold off
plot(t_exact,real(y_exact),'r');
xlabel('time');ylabel('y');axis([t0 tend -inf inf]);
title('Exact solution');
pause
hold on;
plot(t,real(y),'o');
plot(t,real(y));
title('Forward Euler');
