% Declare the relevant parameters
tol=0.01;
h=0.05;
w0=1.0;
tmin=0.0;
tmax=0.5;

nmax=(tmax-tmin)/h;

tlast=tmin;
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)

  t(i) = tlast + h;

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

% Trapezoidal Step
    w_trap_temp = w_euler(i);
    w_trap(i) = wlast_trap + h/2*(deriv(tlast,wlast_trap)+deriv(t(i),w_trap_temp));

  while (abs(w_trap_temp - w_trap(i)) > tol)
      w_trap_temp = w_trap(i);
      w_trap(i) = wlast_trap + h/2*(deriv(tlast,wlast_trap)+deriv(t(i),w_trap_temp));
  end

% Advance time
    wlast_trap =w_trap(i);
    tlast=t(i);

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

end 

% Plot Results for w, y

plot(t,w_trap,'b');
hold on;
plot(t,y,'r');

% Label the plot 

title('Solution to dy/dt=-20y');
xlabel('Time');
ylabel('y');
legend('Trapezoidal','exact')

% 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;
