%% Solving y'=1-t+4*y
T=1;
dt=0.01;


F=@(t,x) 1-t+4*x;

x=1;
time=0;
close all;
FigHandle = figure;
set(FigHandle, 'Position', [100, 100, 1049, 895]);
set(0,'DefaultTextFontSize',30,'DefaultAxesFontSize', 30)

t=0:0.001:T;
plot(t,t/4-3/16+19/16*exp(4*t),'-r'); 
hold all

N=ceil(T/dt);
y=zeros(N+1,1)';
z=y; u=z;
% these are initial conditions at t_0=0
y(1)=1;
z(1)=1;
u(1)=1;

for n=1:N
    % improved Euler
    t=dt*(n-1);
    k1=F(t,z(n));
    k2=F(t+dt, z(n)+dt*k1);
    z(n+1)=z(n)+dt/2*(k1+k2);
    % forward Euler
    y(n+1)=y(n)+dt*F(t, y(n));
    
    
    %% Backward Euler
    u(n+1)=(dt-dt*(t+dt)+u(n))/(1-4*dt); % t is t_n and t+dt is t_{n+1}
end

time=0:dt:dt*N;
hold all;
plot(time, y,'*b','MarkerSize',12);
plot(time, u,'xk','MarkerSize',12);
plot(time, z,'om','MarkerSize',12);
legend('true solution','forward Euler', 'backward Euler','improved Euler')
