% Declare the relevant parameters
h=0.2;
w0=1.0;
tmin=0.0;
tmax=5.0;

nmax=(tmax-tmin)/h;

tlast=tmin;
wlast_euler=w0;
wlast_trap=w0;

t(1) = tmin;
w_euler(1) = w0;
w_trap(1) = w0;
y(1) = w0;

% Beginning of Loop

for i=2:(nmax+1)

% Euler Method
    w_euler(i) = wlast_euler + h*deriv(tlast,wlast_euler);

% Trapezoidal Method
    w_trap(i) = wlast_trap * (1-10*h)/(1+10*h);
%    for dy/dt = + 20y  :    w_trap(i) = wlast_trap * (1+10*h)/(1-10*h);

% Advance time
    t(i) = tlast + h;
    wlast_euler =w_euler(i);
    wlast_trap =w_trap(i);
    tlast=t(i);

% Exact Solution
    y(i) = exact(t(i));

end 

% Exact Solution, fine grid

  tlast=tmin;
  y_fine(1) = w0;
  h=h/10

% Beginning of Loop

  for i=2:(10*nmax+1)
    t_fine(i) = tlast + h;
    y_fine(i) = exact(t_fine(i));
    tlast=t_fine(i);
  end



% Plot Results for w, y

plot(t,w_euler,'b');
%hold on;
%plot(t,w_trap,'og');
%hold on;
%plot(t_fine,y_fine,'r');

% Label the plot 

%title('Solution to dy/dt=t*exp(3*t)-2y');
xlabel('Time');
ylabel('y');
%legend('Euler','Trapezoidal','exact')
legend('Euler')

% Save the results (for each timestep) in a file

file = zeros(size(t,2),3);

file(:,1)=t';
file(:,2)=y';
file(:,3)=w_euler';

save output file -ASCII;
